- © 2020 American Society of Plant Biologists. All rights reserved.
Abstract
Many solanaceous plants secrete acylsugars, which are branched-chain and straight-chain fatty acids esterified to Glu or Suc. These compounds have important roles in plant defense and potential commercial applications. However, several acylsugar metabolic genes remain unidentified, and little is known about regulation of this pathway. Comparative transcriptomics between low- and high-acylsugar-producing accessions of Solanum pennellii revealed that expression levels of known and novel candidate genes (putatively encoding beta-ketoacyl-(acyl-carrier-protein) synthases, peroxisomal acyl-activating enzymes, ATP binding cassette (ABC) transporters, and central carbon metabolic proteins) were positively correlated with acylsugar accumulation, except two genes previously reported to be involved in acylglucose biosynthesis. Genes putatively encoding oxylipin metabolic proteins, subtilisin-like proteases, and other antimicrobial defense proteins were upregulated in low-acylsugar–producing accessions. Transcriptome analysis after biochemical inhibition of biosynthesis of branched-chain amino acids (precursors to branched-chain fatty acids) by imazapyr showed concentration-dependent downregulation of known and most acylsugar candidate genes, but not defense genes. Weighted gene correlation network analysis identified separate coexpressed gene networks for acylsugar metabolism (including six transcription factor genes and flavonoid metabolic genes) and plant defense (including genes putatively encoding NB-ARC and leucine-rich repeat sequences, protein kinases and defense signaling proteins, and previously mentioned defense proteins). Additionally, virus-induced gene silencing of two trichomes preferentially expressed candidate genes for straight-chain fatty acid biosynthesis confirmed their role in acylsugar metabolism.
INTRODUCTION
Plants synthesize a wide variety of specialized metabolites that provide selective advantages in specific environmental conditions. Different classes of specialized metabolites are found in specific taxonomic groups, and many of them are valuable phytochemicals. Acylsugars are nonvolatile and viscous specialized metabolites secreted through trichomes of solanaceous plants (Fobes et al., 1985; Kroumova et al., 2016; Moghe et al., 2017), and their roles in plant defense have been studied extensively. For example, acylsugars from Solanum pennellii act as feeding and oviposition deterrents for insect pests and exert toxic effects on different herbivores (Hawthorne et al., 1992; Juvik et al., 1994; Liedl et al., 1995). Acylsugars from other genera also have important roles in providing protection against herbivores and plant pathogens (Chortyk et al., 1997; Hare, 2005; Luu et al., 2017). Additionally, amphipathic acylsugar molecules are believed to reduce the surface tension of adsorbed dew, thereby providing an additional source of water for the plants (Fobes et al., 1985). Acylsugars also have potential applications as pesticides (Puterka et al., 2003), food additives (Hill and Rhode, 1999), cosmetic and personal care products (Hill and Rhode, 1999), antibiotics (Chortyk et al., 1993), and anti-inflammatory compounds (Pérez-Castorena et al., 2010). Due to these important biological roles and potential commercial applications of acylsugars, a detailed understanding of the genetic network involved in acylsugar metabolism is required for successful crop breeding programs and metabolic engineering of acylsugar production.
In S. pennellii, acylsugars are mostly 2,3,4-tri-O-acylated Glu esters, and some Suc esters, with C4 to C12 fatty acids (Walters and Steffens, 1990; Shapiro et al., 1994; Schilmiller et al., 2015). Predominant branched-chain fatty acids (BCFAs) include 2-methylpropanoic acid, 3-methylbutanoic acid, 2-methylbutanoic acid, and 8-methylnonanoic acid, whereas predominant straight-chain fatty acids (SCFAs) include n-decanoic acid and n-dodecanoic acid. BCFAs are derived from branched-chain amino acids (BCAAs), whereas SCFAs are hypothesized to be derived by a fatty acid synthase (FAS)-mediated de novo fatty acid biosynthetic process (Walters and Steffens, 1990).
Biosynthesis of acylsugars can be divided into two phases: (1) synthesis of fatty acyl chains, and (2) esterification of these acyl molecules to Glu or Suc (Figure 1A). Three acylsugar acyltransferases (ASATs) are involved in phase 2 of acylsucrose biosynthesis (Schilmiller et al., 2015; Fan et al., 2016, 2017). During the preparation of this article, an invertase (acylsucrose fructofuranosidase1; ASFF1) was reported to be capable of producing acylglucose from acylsucrose (Leong et al., 2019). Previous biochemical studies reported a UDP-glucose:fatty acid glucosyltransferase (UDP-Glc:FA GT; Ghangas and Steffens, 1993; Kuai et al., 1997) and a serine carboxypeptidase-like glucose acyltransferase (SCPL GAT; Li et al., 1999; Li and Steffens, 2000) that initiate acylglucose biosynthesis in S. pennellii (Supplemental Figure 1). However, to our knowledge, no genetic evidence is available to support this model.
Current Model of Acylsugar Production in S. pennellii and Workflow to Identify Candidate Genes for Acylsugar Metabolism.
(A) BCFAs are derived from BCAAs, and SCFAs are assumed to be produced by FAS. ASATs use fatty acyl-CoAs for acylsucrose biosynthesis. Recently, an invertase-like enzyme (ASFF1) was reported, which produces acylglucose from acylsucrose. Enzymes are in red font. Unidentified enzymes are indicated with a question mark. Double arrows indicate more than one enzymatic step. C2, 3, 4 indicates esterification at respective positions on Glu. P2, 3, 4 indicates esterification at respective positions on the Suc pyranose ring.
(B) Workflow to identify acylsugar candidate genes. Comparative transcriptomics identified 1,679 DEGs between three low-acylsugar–producing accessions (“LOW”) and three high-acylsugar–producing accessions (“HIGH”). An independent comparison with higher sequencing coverage identified 3,524 DEGs between another low-acylsugar–producing accession, LA1920, and the high-acylsugar–producing accession LA0716. In total, 1,087 DEGs were common to both comparisons. Expression levels of selected candidate genes were determined in LA1302, which accumulates an intermediate level of acylsugars (“MEDIUM”). Analysis of gene expression in response to inhibition of BCAA biosynthesis narrowed the list of candidates, and WGCNA (Langfelder and Horvath, 2008) was performed to identify the coexpressed gene network associated with acylsugar accumulation. Additionally, VIGS led to functional validation of two candidate genes.
Acylsugars are exuded through type-IV trichomes in S. pennellii (Fobes et al., 1985; Slocombe et al., 2008), and many acylsucrose biosynthetic genes are expressed in trichome tip cells (Schilmiller et al., 2012, 2015; Fan et al., 2016). However, analysis of periclinal chimeras indicates that acylglucose biosynthesis in S. pennellii is not trichome-specific (Kuai et al., 1997). Therefore, we used leaf samples for comparative transcriptomics studies reported here. Although this approach underrepresents actual expression differences of trichome-specific genes, it ensures that a complete snapshot of candidate gene networks could be captured.
Different accessions of S. pennellii produce different amounts of total acylsugars (Shapiro et al., 1994). For example, LA0716 produces acylsugars up to 20% of its leaf dry weight (Fobes et al., 1985), whereas in many low-acylsugar–producing accessions, acylsugars make up <1% of leaf dry weight (Shapiro et al., 1994). Here, we exploited this natural genetic variation among S. pennellii accessions to identify candidate biosynthetic and regulatory genes for this trait. Comparative transcriptomics after biochemical inhibition of BCAA biosynthesis further refined the list of candidate genes and illuminated possible regulatory mechanisms of acylsugar production. Functional analysis through virus-induced gene silencing (VIGS) of two candidate genes for SCFA biosynthesis corroborated their involvement in acylsugar metabolism and validated our approach to identifying genes involved in this pathway. Additionally, gene coexpression network analysis revealed genetic networks involved in acylsugar metabolism and plant defense pathways.
RESULTS
Comparative Transcriptomics between Low- and High-Acylsugar–Producing Accessions
We performed differential gene expression analysis between a group of three low-acylsugar–producing accessions (LA1911, LA1912, and LA1926; collectively referred to as “LOW” in this article) and a group of three high-acylsugar–producing accessions (LA1941, LA1946, and LA0716; collectively referred to as “HIGH” in this article) of S. pennellii (Figure 1B; Supplemental Table 1; Shapiro et al., 1994). Three individual plants for each of the six accessions, with an average of >25 million processed paired-end reads for each plant, were used for this analysis (Supplemental Table 2). A total of 19,379 genes passed our filtering criteria for minimum expression levels, and we obtained 1,679 differentially expressed genes (DEGs; see “Methods”; Supplemental Data Set 1: Sheet 1). Of the 1,679 DEGs, 931 were upregulated and 748 were downregulated in the “HIGH” group.
We also compared transcriptomes between another low-acylsugar–producing accession, LA1920, and the high-acylsugar–producing accession LA0716, with each accession having four biological replicates and slightly higher sequencing coverage, as an independent examination of differential gene expression. In this analysis, 19,967 genes passed our filtering criteria for minimum expression levels (see “Methods”), and we identified 3,524 DEGs (1,465 upregulated and 2,059 downregulated genes in LA0716; Supplemental Data Set 1: Sheet 2). To identify candidate genes most likely to be involved in acylsugar metabolism, we compared “LOW” versus “HIGH” DEGs with LA1920 versus LA0716 DEGs (Figure 1B). A total of 1,087 DEGs were common to both comparisons (586 upregulated and 501 downregulated genes in high-acylsugar–producing accessions; Supplemental Data Set 1: Sheet 3). Enrichment analysis revealed significant overrepresentation of gene ontology (GO) terms such as “secondary metabolite biosynthetic process” (GO:0044550), “fatty acid biosynthetic process” (GO:0006633), “branched-chain amino acid biosynthetic process” (GO:0009082), and “defense response” (GO:0006952) in the list of common DEGs (Supplemental Figure 2).
We also examined the leaf transcriptome of accession LA1302 (three biological replicates), which accumulates an intermediate level of acylsugars (referred to as “MEDIUM” in this article; Supplemental Table 1; Shapiro et al., 1994), to determine whether expression levels of candidate genes were consistent with the amount of acylsugar production (Figure 1B). Multidimensional scaling plots showed a clear distinction in gene expression profiles among groups of accessions that accumulate different amounts of acylsugars (Supplemental Figures 3A and 3B). Volcano plots showed high log2(fold-change) values associated with low false discovery rate (FDR) values, indicating robust differential expression (Supplemental Figures 3C and 3D). On the other hand, 10 “housekeeping genes” showed similar expression levels across all biological groups (Supplemental Figure 4).
BCAA Metabolic Genes Are Upregulated in High-Acylsugar–Producing Accessions
Because BCAAs are precursors to acylsugar BCFAs (Walters and Steffens, 1990), we examined genes involved in BCAA metabolism. Many genes involved in both biosynthesis of BCAAs in plastids and conversion of BCAAs to branched-chain acyl molecules in mitochondria (Binder, 2010; Maloney et al., 2010) were upregulated in high-acylsugar–producing accessions (Table 1; Supplemental Figure 5). Most of these BCAA/BCFA metabolic DEGs showed intermediate levels of expression in the “MEDIUM” accession LA1302 (Figure 2A). Correlation analysis (Spearman’s rank correlation coefficient, SRCC; denoted by ρ) using gene expression profiles in 29 samples revealed a strong positive correlation among most of these genes (Figure 3).
Heatmaps Showing Expression Levels of Genes with Known and Putative Functions in Acylsugar Metabolism. Genes are designated by their gene identifier numbers (Sopen IDs).
(A) Acylsugar phase-1–related genes (Table 1). Expression levels of many candidate genes in LA1302 (“MEDIUM” accession) were intermediate between low- and high-acylsugar–producing accessions.
(B) Acylsugar phase-2–related genes (Table 2). Sopen04g001210 (marked with an asterisk) was identified as a carboxylesterase gene related to ASH (Schilmiller et al., 2016). Sopen03g040490 (marked with a double-asterisk) has been recently reported as a trichome-expressed invertase gene that is capable of producing acylglucose from acylsucrose (Leong et al., 2019). AG, acylglucose phase 2 (previous model).
(C) Other genes related to acylsugar metabolism (Table 3). Two genes (Sopen04g023150 and Sopen01g047950) were removed from the heatmap due to their very high expression levels. Sopen04g023150 had FPKM values of 570, 624, 455, 179, and 254 in “LOW,” LA1920, LA1302, LA0716, and “HIGH” groups, respectively. Sopen01g047950 had FPKM values of 61, 203, 79, 29, and 19 in “LOW,” LA1920, LA1302, LA0716, and “HIGH” groups, respectively. AAE1, acyl-activating enzyme 1; ABC, ABC transporters; AG, acylglucose phase 2 (previous model); ASAT, acylsucrose acyltransferase; ASH, acylsugar hydrolase; BCAA, branched-chain amino acid; CCM, central carbon metabolism; FAS, fatty acid synthase components; FPKM, fragments per kilobase of transcript per million mapped reads.
Correlation Among Expression Profiles of Selected Genes.
FPKM values for 29 phase-1- and phase-2–related genes (from Tables 1 and 2, respectively) in our 29 samples were used to determine pairwise SRCC. A [29 × 29] matrix was created to visualize correlation. Red, white, and blue colors indicate positive, zero, and negative correlations, respectively. Darker colors indicate stronger correlations. SRCC values and associated P values are given in Supplemental Data Set 2. Sopen04g001210 (marked with an asterisk) is a carboxylesterase gene related to ASH. Sopen03g040490 (marked with a double-asterisk) is the recently identified invertase gene. Abbreviations: AG, acylglucose phase 2 (previous model).
Genes Putatively Encoding a KAS IV/KAS II-like Enzyme, Other FAS Components, and AAEs Are Upregulated in High-Acylsugar–Producing Accessions
Sopen12g004230, predicted to encode a beta-ketoacyl-(acyl-carrier-protein) synthase II (KAS II)-like enzyme, had 39-fold higher expression in the “HIGH” group (FDR = 4.40E−18). The SCFA profile of S. pennellii acylsugars (predominantly C10 and C12, with a trace amount of C8 in LA1302) is strikingly similar to that of the seed oil in many Cuphea species (Shapiro et al., 1994; Dehesh et al., 1998; Slabaugh et al., 1998; Schütt et al., 2002). A specialized version of KAS II, referred to as KAS IV, is an important determinant of chain length in Cuphea medium-chain (C8 to C12) fatty acids (Dehesh et al., 1998; Slabaugh et al., 1998; Schütt et al., 2002). We compared amino acid sequences of five KAS IV/KAS II-like enzymes from four Cuphea species with the Sopen12g004230 sequence, but no similarity was found. However, sequence analysis of the adjacent Sopen12g004240 coding region (1,785-amino acids; annotated as reverse transcriptase) revealed a KAS domain at the C-terminal end (1,374 to 1,784), and BLAST (https://blast.ncbi.nlm.nih.gov/Blast.cgi) analysis of Cuphea KAS IV/KAS II-like enzymes showed high sequence similarity with Sopen12g004240 (77 to 83% identity, 87 to 92% similarity, 82 to 90% coverage, e-value 0.0; Supplemental Table 3). Apparently, a transposon insertion into an intron of Sopen12g004230-Sopen12g004240 led to misannotation of this locus as two separate genes (Supplemental Figure 6A). We confirmed that this locus produces a single transcript with reverse transcription (RT)-PCR analysis, and Sanger sequencing confirmed that this single transcript can produce a functional KAS protein (Supplemental Figures 6B and 6C). An intraspecific F2 population derived from the cross between accessions LA0716 and LA1912 identified a quantitative trait locus on the upper arm of chromosome 12 that is associated with C10 SCFA production (Blauth et al., 1999). The locus containing Sopen12g004230-Sopen12g004240 is in this region.
DEGs putatively encoding a KAS III, other components of the FAS complex, and acyl-activating enzymes (AAEs; ATP/AMP binding proteins that form fatty acyl-CoA molecules from free fatty acids, ATP and CoA; Shockey et al., 2003) were also upregulated in high-acylsugar–producing accessions (Table 1). Expression profiles of these genes showed a very strong positive correlation among themselves, as well as with BCAA/BCFA metabolic DEGs (Figure 3; Supplemental Data Set 2).
Acylsugar Phase-2 Metabolic Genes Are Upregulated in High-Acylsugar–Producing Accessions and Show a Strong Positive Correlation with Phase-1 DEGs
Genes encoding ASATs and the recently identified invertase (SpASFF1), which are involved in phase 2 of acylsugar biosynthesis (Schilmiller et al., 2015; Fan et al., 2016, 2017; Leong et al., 2019), were upregulated in high-acylsugar–producing accessions (Table 2). ASATs use acyl-CoA molecules as their substrates (Schilmiller et al., 2012, 2015; Fan et al., 2016, 2017), and free SCFAs produced by FAS-catalyzed de novo fatty acid biosynthesis must be activated to their acyl-CoA derivatives (SCFA-CoA) before they can be used by ASATs (Figure 1A). All three ASAT genes exhibited very strong positive correlation with both acyl-activating enzyme1 (AAE1) genes and FAS components (Figure 3; Supplemental Data Set 2).
Recently, three genes encoding acylsugar hydrolases (ASHs; carboxylesterases that remove acyl groups from acylsucrose molecules) and a related carboxylesterase gene (Sopen04g001210) were reported in cultivated and wild tomato (Schilmiller et al., 2016). Expression of these genes, except Sopen05g030120 (Sp-ASH1) showed strong positive correlation with genes putatively involved in phase-1 acyl chain synthesis (Figure 3; Supplemental Data Set 2). The expression profile of ASH1 in different tissues of Solanum lycopersicum led the authors to propose additional non-trichome–localized functions for Sl-ASH1 (Schilmiller et al., 2016); our results further support the idea that ASH1 may not have a critical role in trichome acylsucrose metabolism.
Expression Levels of Two Previously Reported Acylglucose Biosynthetic Genes Are Not Consistent with Acylsugar Levels
Expression profiles of two genes encoding UDP-Glc:FA GT and SCPL GAT, respectively (previous model of acylglucose biosynthesis; Supplemental Figure 1; Ghangas and Steffens, 1993; Kuai et al., 1997; Li et al., 1999; Li and Steffens, 2000) were not positively correlated with acylsugar levels (Table 2; Figure 2) and other acylsugar metabolic genes (Figure 3; Supplemental Data Set 2).
Genes Putatively Encoding ATP Binding Cassette Transporters and Central Carbon Metabolic Proteins Are Upregulated in High-Acylsugar–Producing Accessions
Five genes belonging to the ATP-binding cassette (ABC) transporter family showed differential expression, and three of them (Sopen12g034820, Sopen04g005380, and Sopen03g001870) were upregulated in high-acylsugar–producing accessions (Table 3). Recently, the role of central carbon metabolism (CCM) in supporting specialized metabolism in glandular trichomes of S.lycopersicum and Solanum habrochaites was investigated, and it was concluded that glandular trichomes import most of their fixed carbon from underlying leaf tissues, despite being able to perform some photosynthetic activities (Balcke et al., 2017). Special roles of chloroplast malic enzyme in supplying pyruvate and glutathione s-transferase (GST) in alleviating oxidative stress were also reported. We did find genes putatively encoding these central carbon metabolic (CCM) proteins in the list of 1,087 DEGs (Table 3; Figure 2C).
Trichome-Enriched Expression of Acylsugar Candidate Genes
Genes encoding ASATs and the recently identified invertase (SpASFF1) are expressed in apical cells of type-I/IV glandular trichomes (Schilmiller et al., 2012, 2015; Fan et al., 2016; Leong et al., 2019). Using data of Ning et al. (2015), who published expression profiles of S. lycopersicum genes in isolated stem trichomes versus shaved stems, we investigated if S. lycopersicum homologs of candidate genes show a trichome-enriched expression pattern. We used reciprocal BLAST to identify putative S. lycopersicum orthologs of candidate genes (Supplemental Data Set 3), and they showed trichome-enriched expression pattern (26- to 562-fold for AAE1 members, 143- to 470-fold for ABC transporter genes, and 16- to 287-fold for CCM genes; Supplemental Table 4).
Genes Putatively Encoding Oxylipin Metabolic Proteins and Other Defense Proteins Are Upregulated in Low-Acylsugar–Producing Accessions
Defense response genes were present in the list of 1,087 DEGs (GO:0006952; Supplemental Figure 2), and we found that genes putatively involved in oxylipin biosynthesis (fatty acid desaturase [FAD] and lipoxygenase [LOX]), metabolism (oxo-phytodienoate reductase and GST), and signaling pathways (TGA transcription factors; Mueller et al., 2008), as well as those putatively encoding subtilisin-like proteases (SLPs) and additional defense proteins were upregulated in low-acylsugar–producing accessions (Table 4; Figure 4). Based on expression profiles, these putative defense genes showed moderate to strong correlation with each other (Supplemental Figure 7), indicating an active network of innate immunity in low-acylsugar–producing accessions. According to data of Ning et al. (2015), putative S. lycopersicum orthologs of these defense DEGs do not show trichome-enriched expression (except two of 28 DEGs; Sopen09g006360 and Sopen10g019020; Supplemental Table 4).
Heatmap Showing Expression Levels of Genes Putatively Involved in Plant Defense (Table 4 Genes).
Expression levels of many genes in LA1302 were similar to low-acylsugar–producing accessions but different from high-acylsugar–producing accessions. Some oxylipin-related defense genes, which were upregulated in low-acylsugar–producing accessions, had higher expression levels than amine oxidase (AOX)-related defense genes, which were upregulated in high-acylsugar–producing accessions. OLM, oxylipin metabolism.
Amine Oxidase-Related Defense Response Genes Are Upregulated in High-Acylsugar–Producing Accessions
Genes putatively encoding amine oxidase domains (FLOWERING LOCUS D-like proteins; required for systemic acquired resistance in Arabidopsis [Arabidopsis thaliana]; Singh et al., 2013) and a corresponding aldehyde oxidase (Cona et al., 2006) were upregulated in high-acylsugar–producing accessions (Table 4). H2O2 generated by amine oxidase- and aldehyde oxidase-catalyzed reactions acts as a signal molecule and plays an important role during pathogen invasion by acting on extensin proteins (Niebel et al., 1993; Jackson et al., 2001; Cona et al., 2006). Two genes predicted to encode extensin-like glycoproteins were also upregulated in high-acylsugar–producing accessions (Table 4), and one of them (Sopen05g003410) showed strong positive correlation with amine oxidase- and aldehyde oxidase-encoding genes (Supplemental Figure 7). Although amine oxidase-related defense genes were upregulated in high-acylsugar–producing accessions, it is important to note that their expression levels were much lower than some oxylipin-related defense genes in low-acylsugar–producing accessions (Figure 4). Among these amine oxidase-related DEGs, we identified a putative S. lycopersicum ortholog (Solyc01g005850) for only one gene (Sopen01g001850; extension-like protein), and according to the data of Ning et al. (2015), Solyc01g005850 does not show trichome-enriched expression (Supplemental Table 4).
Biochemical Inhibition of BCAA Biosynthesis Changes Transcript Levels of Many Acylsugar Candidate DEGs, But Not Defense DEGs
Acetolactate synthase (ALS) catalyzes the first common step in the biosynthesis of BCAAs (Supplemental Figure 5A), and biochemical inhibition of this enzyme significantly lowers acylsugar production in S. pennellii (Walters and Steffens, 1990). To determine whether inhibition of this enzyme changed expression of ALS and other candidate genes involved in acylsugar metabolism, we treated leaves of S. pennellii (the high-acylsugar–producing accession LA0716) with the ALS inhibitor imazapyr at 0.1 mM and 1 mM concentrations.
Expression of 171 of the 1,087 DEGs previously identified by comparative transcriptomics was affected by imazapyr treatment at both concentrations, including many of the DEGs involved in BCAA/BCFA metabolism (Table 5; Supplemental Data Set 4: Sheet 1). It should be noted that at 1 mM imazapyr, expression of Sopen04g026270 and Sopen01g028100 (branched-chain keto acid dehydrogenase complex E1 and E2 subunits, respectively) was significantly increased, possibly because the amount of 2-ketobutyrate, a substrate of ALS that also can be used directly by the dehydrogenase complex to make propionyl-CoA, builds up (Walters and Steffens, 1990).
DEGs putatively involved in SCFA metabolism (FAS componentscand AAE1 members) also showed significant downregulation at both 0.1 mM and 1 mM imazypyr in a concentration-dependent manner. AAE1 proteins have a peroxisomal location in Arabidopsis (Reumann et al., 2009); one gene (Sopen11g007710) predicted to encode a PMP22/Mpv17 family peroxisomal membrane protein had 7.7-fold (FDR = 1.17e−08) and 165-fold (FDR = 3.62E−24) lower expression at 0.1 mM and 1 mM imazapyr treatment, respectively (Supplemental Data Set 4: Sheet 1). Sopen11g007710 was upregulated 29-fold in the “HIGH” group (FDR = 2.80E–11; Supplemental Data Set 1), and its putative ortholog in S. lycopersicum (Solyc11g012990; Supplemental Data Set 3) has 1,420-fold higher expression in isolated trichomes than shaved stems, consistent with trichome-enriched expression of AAE1 members (Ning et al., 2015). PMP22 is presumably involved in controlling permeability of peroxisomal membrane (Brosius et al., 2002). These results suggest a role of peroxisome in acylsugar metabolism.
DEGs encoding ASATs and the invertase showed significant decrease in gene expression level in a concentration-dependent manner (by as much as 94-fold, FDR = 1.65E–14), as did ASHs and the related carboxylesterase Sopen04g001210 (by as much as 70-fold, FDR = 8.82E–50; Table 5). On the other hand, Sopen01g049990 and Sopen10g020280 (encoding UDP-Glc:FA GT and SCPL GAT, respectively; previous model of acylglucose biosynthesis) showed slightly higher and similar expression levels, respectively, in response to imazapyr treatment.
DEGs putatively encoding CCM proteins and three ABC transporters were significantly downregulated (by as much as 45-fold, FDR = 4.99E−10 for RUBISCO small subunit; Table 5). Imazapyr treatment in Arabidopsis resulted in significant induction of nine genes (out of the 10 genes reported) encoding ABC transporters, indicating their role in detoxification process (Manabe et al., 2007). Similarly, 15 putative ABC transporter genes in S. pennellii were upregulated in response to imazapyr either at higher concentration or at both concentrations (Supplemental Table 5). However, three putative ABC transporter DEGs were repressed by imazapyr treatment, consistent with expression profiles of other acylsugar metabolic genes.
Unlike known and candidate acylsugar metabolic DEGs, amine oxidase-related defense DEGs did not show concentration-dependent downregulation in response to ALS inhibition (Supplemental Data Set 4: Sheet 2). The overall pattern of gene expression in response to imazapyr suggests that additional genes involved in acylsugar metabolism may be among these 171 DEGs. Of the 13 transcription factor DEGs in the initial list of 1,087, nine were present in the imazapyr treatment list (Table 3; Figure 2C). Putative S. lycopersicum orthologs of six of these nine DEGs have a trichome-enriched expression profile (Supplemental Table 4; Ning et al., 2015).
Identification of Coexpressed Gene Network Associated with Acylsugar Metabolism
Using expression profiles of S. pennellii genes in our 38 samples (29 from different accessions and nine from the imazapyr treatments), weighted gene correlation network analysis (WGCNA; Langfelder and Horvath, 2008) identified 41 coexpressed gene modules. Known and candidate acylsugar metabolic genes were clustered in the “darkorange” module (182 genes; Figure 5A; Supplemental Figure 8). One-third of the genes that responded to imazapyr (57/171) were also placed in this module (Supplemental Data Set 4: Sheet 3). We selected the top-50 most strongly connected genes (based on intramodular connectivity as determined by WGCNA) for each of the three ASAT genes and the invertase gene, and merged them to identify the most strongly interconnected gene network. This network included genes involved in BCAA/BCFA metabolism, SCFA synthesis (FAS components) and their activation (AAEs and the related peroxisomal membrane protein PMP22/Mpv17), acylsucrose metabolism (three ASATs, ASH2, ASH3 and the related carboxylesterase Sopen04g001210), as well as genes putatively encoding three ABC transporters, CCM proteins, six transcription factors (including three AP2-family transcription factors), and other proteins (Figure 5B).
WGCNA (Langfelder and Horvath, 2008).
(A) Association of module eigengenes (MEs) with expression profiles of nine selected genes. WGCNA identified 41 modules; for each module, an ME was selected as representative of the expression profile for that module (first principle component). Numbers indicate SRCC values between MEs of two selected modules and nine genes’ expression profiles. P values associated with SRCC are given in parentheses (9E-23 = 9 × 10−23). Blue, white, and red colors indicate negative, zero, and positive correlations, respectively. Acylsugar metabolic genes (ASATs and invertase) and LOX-correlated genes were clustered in “darkorange” and “green” modules, respectively. NB-ARC-LRR indicates a recently characterized sequence (Sopen09g036080) that confers resistance against tospoviruses in S. lycopersicum and wild tomato species (Zhu et al., 2017). A complete list of WGCNA modules is given in Supplemental Figure 8.
(B) Simplified acylsugar metabolic gene network. Nodes and edges represent genes and intramodular connectivities, respectively. We used WGCNA to identify the top-50 genes most strongly connected to each of the three ASAT genes and the invertase gene. These four top-50 lists contained only 58 different genes.
However, not all genes in this network had an obvious connection to acylsugar metabolism. One gene on chromosome 11 (Sopen11g003320; UDP-glucose:catechin glucosyltransferase; Noguchi et al., 2008) and three sequential genes on chromosome 6 (Sopen06g034810, Sopen06g034820, and Sopen06g034830; myricetin methyltransferase (Schmidt et al., 2012; Kim et al., 2014) showed strong intramodular connectivity with acylsugar metabolic genes. These flavonoid metabolic genes were also strongly downregulated in response to imazapyr treatment (Supplemental Figure 9), and their putative orthologs in S. lycopersicum are preferentially expressed in trichomes (Ning et al., 2015). In addition to abundant acylsugar compounds, flavonoid compounds, such as methylated myricetin, have been reported in S. pennellii trichomes (McDowell et al., 2011). Together, these results suggest that metabolisms of acylsugar and flavonoid compounds are strongly connected.
Plant Defense Gene Network
We selected Sopen01g002520 (LOX) as a representative of the oxylipin-mediated defense system, due to its robust differential expression (Figure 4; Table 4), and WGCNA identified coexpressed genes in the “green” module. Enriched GO terms in this module include “defense response” (GO:0006952), “oxylipin biosynthetic process” (GO:0031408), and “glutathione metabolic process” (GO:0006749; Figure 6; Supplemental Figure 10). In the plant innate immune system, recognition and successful prevention of pathogen invasion rely on a series of molecular events, including pathogen detection at the cell surface by receptor kinases, a series of phosphorylation events with cytoplasmic kinases, nucleotide binding, release of intracellular calcium, protein ubiquitination, and transcriptional reprogramming (Couto and Zipfel, 2016). GO terms associated with these processes were significantly enriched in the “green” module (FDR = 8.83E−17 to 4.6E−02). Twenty of 104 expressed S. pennellii NB-ARC sequences and 23 of 150 leucine-rich repeat (LRR) sequences were clustered in the “green” module (Supplemental Data Set 5), including the recently characterized NB-ARC-LRR sequence Sw-5b (Sopen09g036080) that confers broad-spectrum resistance against American-type tospoviruses in S. lycopersicum and wild tomato species (Zhu et al., 2017).
Enrichment of Selected GO Terms Associated with WGCNA “Green” Module Genes.
Test set and reference set indicate “green” module genes and S. pennellii annotated sequences, respectively. Enriched sequences include those putatively involved in oxylipin metabolism and plant defense signaling (Couto and Zipfel, 2016). Enrichment analysis was performed using Blast2GO software (Fisher’s exact test). Only significant GO terms are shown (FDR < 0.05). A complete list of all significant GO terms, GO IDs, and associated FDR and P values is given in Supplemental Figure 10.
Functional Validation of Two Candidate Genes Involved in SCFA Biosynthesis
To our knowledge, there are no reports of genes involved in acylsugar SCFA synthesis. We selected two candidate genes for functional validation (Sopen08g002520 and Sopen12g004240; predicted to encode KAS III and KAS IV/KAS II-like enzymes, respectively). These genes have 235- and 60-fold, respectively, higher expression levels in isolated stem trichomes compared with underlying tissue (Figure 7A). In fact, Sopen08g002520 showed a SpASAT1-like trichome-enriched expression profile. We identified Solyc08g006560 and Solyc12g009260 as putative S. lycopersicum orthologs of Sopen08g002520 and Sopen12g004240, respectively, and these two tomato orthologs also show trichome-enriched expression profiles (90- and 103-fold, respectively; Supplemental Table 4; Ning et al. (2015)). Phylogenetic analyses supported these orthologous relationships and revealed distinct clades of KAS III and KAS II sequences in solanaceous species (Supplemental Figure 11; Supplemental Files 1 to 4).
Functional Validation of Two Candidate KAS Genes.
(A) Transcript levels of Sopen08g002520 (KAS III) and Sopen12g004240 (KAS IV/II-like) in isolated stem trichomes and underlying tissues (shaved stems) were measured by RT-qPCR. Error bars indicate se (n = 5 individual plants).
(B) Acyl chain composition of acylsugars extracted from leaves of control plants (empty TRV vector) and two other groups of plants in which Sopen08g002520 and Sopen12g004240, respectively, were targeted for VIGS. VIGS of both Sopen08g002520 and Sopen12g004240 resulted in significant decrease in SCFAs (n-decanoic acid and n-dodecanoic acid). Relative abundance of acyl chains was measured by GC–MS. Asterisks denote statistical significance compared with control plants as determined by one-way ANOVA (Dunnett’s test; *P < 0.05, **P < 0.01, ***P < 0.001). Error bars indicate se (n = 10 individual plants).
We targeted these two S. pennellii loci in LA0716 for VIGS using tobacco rattle virus (TRV)-based silencing vectors, and acylsugar SCFA production was reduced by up to 40% in both n-C10 and n-C12 amount (P < 0.001; Figure 7B). VIGS resulted in significant downregulation of target genes (60% and 67% reduction in transcript levels for Sopen08g002520 and Sopen12g004240, respectively; Supplemental Figures 12A and 12B). Residual expression presumably is due to the fact that VIGS in S. pennellii is incomplete and inconsistent (Supplemental Figure 12C). Additionally, no morphologically distinct areas would allow for better area-selective metabolite profiling, because no significant morphological differences were observed between silenced and control plants (Supplemental Figure 12D). Nevertheless, our results indicate that the two candidate genes, which are preferentially expressed in trichomes, play important roles in acylsugar SCFA biosynthesis. This is consistent with the previous prediction that KAS enzyme(s) other than KAS I might be involved in SCFA production in S. pennellii acylsugars (Slocombe et al., 2008).
DISCUSSION
Acylsugars are powerful natural pesticides (Puterka et al., 2003), and increasing acylsugar-mediated insect pest resistance has been a target of various breeding programs in solanaceous crops (Bonierbale et al., 1994; Lawson et al., 1997; Alba et al., 2009; Leckie et al., 2012). The success of these breeding programs and bioengineering of acylsugar production will largely depend on unraveling the network of structural and regulatory genes. Here, we used transcriptomics approaches to identify candidate acylsugar and defense gene networks in S. pennellii.
Defense Pathway Genes
Enriched GO terms associated with pathogen recognition and downstream signaling cascades in the “green” module indicate an active defense network in low-acylsugar–producing accessions (Figure 6). Oxylipins (Halitschke and Baldwin, 2003; Prost et al., 2005) and SLPs (Figueiredo et al., 2014) have important roles in plant defense, and higher expression levels of these defense genes predominantly in low-acylsugar–producing accessions (Table 4) may indicate a compensatory mechanism for diminished defense activities of acylsugars in these accessions. Alternatively, expression of these defense genes may be suppressed in high-acylsugar–producing accessions due to the protection provided by acylsugars. Our data do not allow us to distinguish between these two possibilities, but it is important to note that all plants used in this study were grown in the same growth chamber, and no insects or symptoms of disease were apparent on any plant.
Most defense DEGs (both oxylipin-related and amine oxidase-related) had similar expression levels in the “MEDIUM” accession LA1302 and in low-acylsugar–producing accessions (Figure 4). Multidimensional scaling plots also showed that the overall gene expression profile of LA1302 is similar to that of low-acylsugar–producing accessions, but very different from that of high-acylsugar–producing accessions (Supplemental Figures 3A and 3B). These results are consistent with the prediction that LA1302 is more closely related to low-acylsugar–producing accessions (Shapiro et al., 1994). Despite the similarity of defense DEGs expression and the overall transcriptome profile between LA1302 and low-acylsugar–producing accessions, many DEGs with known and putative functions in acylsugar metabolism showed intermediate expression levels in LA1302 (Figure 2). Together, these results suggest that there is no strict linear relationship between acylsugar accumulation and expression of defense pathways.
Acylsugar Metabolic Genes
Our initial hypothesis for this work was that expression levels of previously reported acylsugar biosynthetic genes would be consistent with acylsugar levels in different accessions of S. pennellii, and that additional, unknown genes involved in acylsugar metabolism would be coordinately expressed with them. Our approach was validated by the identification of several DEGs with known functions in acylsugar biosynthesis, such as BCAA/BCFA metabolic genes and ASATs. Next, we used a chemical genetics approach to refine the list of acylsugar candidate genes by treating S. pennellii leaves with imazapyr to inhibit BCAA biosynthesis. This approach allowed us to separate the acylsugar metabolic DEGs from defense DEGs, which did not respond to imazapyr treatment (Table 5; Supplemental Data Set 4). Although we expected that genes related to synthesis of acylsugar BCFAs (derivatives of BCAAs) would be downregulated as a result of lack of substrate availability in response to imazapyr, our results showed that genes related to SCFA metabolism and other aspects of acylsugar metabolism (including known acylsucrose metabolic genes and candidate genes) also responded to imazapyr in a dose-dependent manner. It is important to note that the repression of acylsugar metabolic genes we observed after imazapyr treatment is distinct from the induction of general stress response genes that can result from amino acid starvation (Zhao et al., 1998).
Strong correlations among expression profiles of selected acylsugar metabolic genes (Figure 3) encouraged us to perform gene coexpression network analysis, which identified candidate gene networks involved in acylsugar metabolism and also plant defense (Figures 5 and 6; Supplemental Data Set 5). Additionally, data of Ning et al. (2015) indicate that putative S. lycopersicum orthologs of most of our acylsugar candidate genes, but not defense DEGs, are preferentially expressed in trichomes (Supplemental Table 4). Together, these results strongly support the candidacy of novel acylsugar metabolic genes that we have identified; candidate genes include those encoding KAS enzymes and other FAS components (for SCFA synthesis), AAE1 members (for activating SCFAs for ASATs), a PMP22/Mpv17 family peroxisomal membrane protein (for controlling permeability of peroxisomal membrane), ABC transporters (for acylsugar export), and CCM proteins (for providing carbon to support acylsugar production), including a glutathione s-transferase (for alleviating oxidative stress).
In plant primary metabolism, the condensation reactions of fatty acid biosynthesis are catalyzed by KAS enzymes of the FAS complex, and at least three different KAS isoenzymes (KAS I, KAS II, and KAS III) are involved in de novo fatty acid biosynthesis. We selected two candidate KAS genes that are preferentially, if not exclusively, expressed in S. pennellii trichomes, and knockdown of these two genes led to significant decrease in SCFA synthesis (Figure 7), demonstrating their role in acylsugar metabolism. Duplication of a primary metabolic gene followed by neofunctionalization (Ning et al., 2015) and spatial gene regulation have important functions in regulating trichome metabolism of acylsugar molecules, and KAS sequences also appear to follow that trend.
Regulation of Acylsugar Production
Both comparative transcriptomics and imazapyr treatment studies revealed that transcript levels for the two previously reported phase-2 acylglucose biosynthetic enzymes (UDP-Glc:FA GT and SCPL GAT) were not positively correlated with acylsugar levels (Tables 2 and 5). Recently, Leong et al. (2019) clearly demonstrated the role of the trichome-expressed invertase (SpASFF1) in acylglucose biosynthesis in S. pennellii trichomes. Their results contradict the involvement of UDP-Glc:FA GT and SCPL GAT in the biosynthesis of acylglucose molecules, and our results are consistent with this finding. Our results also suggest that a majority of the acylsugar metabolic pathway, from acyl chain synthesis to acylsugar secretion, represents a remarkable coexpressed genetic network. WGCNA revealed that six transcription factor genes form a strong gene coexpression network with many acylsugar and flavonoid metabolic genes. Downregulation of these transcription factor genes in response to imazapyr treatment parallels significant downregulation of most acylsugar and flavonoid metabolic genes (Table 5; Supplemental Data Set 4). Additionally, putative orthologs of these transcription factor genes in S. lycopersicum exhibit trichome-enriched expression, similar to many genes involved in acylsugar and flavonoid metabolism (Ning et al., 2015). Together, these results suggest important roles for these transcription factors in specialized metabolic pathways of S. pennellii trichomes.
METHODS
Plant Materials
Seeds of all Solanum pennellii accessions were obtained from the C.M. Rick Tomato Genetics Resource Center (University of California, Davis). Seeds were treated with a solution of 20% (v/v) commercial bleach (1.2% [w/v] sodium hypochlorite) for 20 min, and then germinated on moistened filter paper. Seedlings were transferred to soil at the cotyledon stage, and plants were grown in a growth chamber (24°C day/20°C night temperature, 150 mMol m−2 s−1 photosynthetically active radiation, 16-h photoperiod, and 75% humidity).
RNA Sequencing Library Preparation and Sequencing
Three biological replicates were used for each accession in the “LOW” versus “HIGH” comparison, and for the “MEDIUM” accession LA1302. Four biological replicates were used for each accession in the LA1920-versus-LA0716 comparison. Before RNA extraction, secreted acylsugars were removed from the leaf surface of 10-week–old plants by dipping them in ethanol for 2 to 3 s. Leaves were then immediately frozen in liquid nitrogen and stored at −80°C until further use. Total RNA was extracted from leaves using the RNAqueous Total RNA Isolation Kit (Thermo Fisher Scientific). Genomic DNA was removed using the TURBO DNA-Free Kit (Thermo Fisher Scientific). Quality of each RNA sample was analyzed with the Agilent 2200 TapeStation software A01.04 (Agilent Technologies). RNA sequencing (RNA-Seq) libraries were prepared from polyA+-selected RNA samples using the TruSeq RNA Library Preparation Kit v2 (Illumina), and then sequenced on the HiSeq 2500 v4 125-bp × 125-bp paired-end sequencing platform (High Output Mode; Illumina) following the manufacturer’s specifications. Sequencing was performed at the Texas A&M Genomics and Bioinformatics Service Center. Sequence cluster identification, quality prefiltering, base calling, and uncertainty assessment were done in real time using Illumina’s HCS v2.2.68 and RTA v1.18.66.3 software with default parameter setting. Sequencer.bcl basecall files were demultiplexed and formatted into .fastq files using bcl2fastq v2.17.1.14 script conFigureBclToFastq.pl. Sequencing reads were submitted to the National Center for Biotechnology Information Sequence Read Archive under accession number SRP136022.
Quality Control and Mapping of Sequencing Reads
Phred quality score distributions of sequencing reads were analyzed with the software FastQC v0.11.4 (Andrews, 2010). All RNA-Seq libraries had average Phred quality scores of >35, indicating >99.97% base call accuracy. To remove all the low-quality, ambiguous “N” nucleotides (Phred quality score = 2, ASCII code = #), the following settings were applied for quality control of RNA-Seq reads using the program Trimmomatic v0.35 (Bolger et al., 2014b): ILLUMINACLIP:TruSeq3-PE-2.fa:2:30:10, SLIDINGWINDOW = 1:14, MINLEN = 80. Using these settings, >98% of the sequencing reads could be retained, and it ensured that none of the nucleotides in the reads were ambiguous.
Processed reads were mapped to the S. pennellii genome v2.0 (LA0716; Bolger et al., 2014a) using TopHat2 v2.1.0 (Kim et al., 2013). The following parameters were used for the mapping process: -mate-inner-dist = 0, -mate-std-dev = 50, -read-realign-edit-dist = 0, -read-edit-dist = 4, -library-type = fr-unstranded, -read-mismatches = 4, -bowtie-n, -min-anchor-len = 8, -splice-mismatches = 0, -min-intron-length = 50, -max-intron-length = 50,000, -max-insertion-length = 3, -max-deletion-length = 3, -max-multihits = 20, -min-segment-intron = 50, -max-segment-intron = 50,000, -segment-mismatches = 2, -segment-length = 25. A summary of TopHat2 alignment results are given in Supplemental Table 2. Reads mapped to selected loci of interest were visualized with IGV (Thorvaldsdóttir et al., 2013).
DEG Analyses
Uniquely-mapped reads from TopHat2 were counted for all S. pennellii annotated genes using HTSeq-Count (Anders et al., 2015). The following parameters were used for the counting process: -f bam, -r name, -s no, -m union, -a 20. These count files were used for identification of DEGs using the program edgeR (Robinson et al., 2010). Genes with very low expression levels were filtered out before differential testing. For comparison between the “LOW” and “HIGH” groups, genes with >1 count per million (CPM) in at least six samples were retained for analysis after read counts were adjusted for transcript lengths. Library sizes were recalculated after the filtering process, and all the samples had post-filter library sizes of >99.70%. Tagwise dispersions were calculated (prior.df = 3), and genes were identified as differentially expressed when P values, corrected for multiple testing using Benjamini–Hochberg multiple testing correction (Benjamini and Hochberg, 1995), were <0.05 (FDR < 0.05), and had at least 2-fold expression differences. For comparison between LA1920 and LA0716, we used the same pipeline, except genes with >1 CPM in at least two samples were retained for analysis, and tagwise dispersions were calculated with prior.df = 9. In this analysis, all the samples showed >99.85% post-filter library sizes. Lists of DEGs as determined by edgeR (Robinson et al., 2010) are given in Supplemental Data Set 1.
Functional Annotation of DEGs
Transcript sequences of DEGs were extracted from the annotated gene models, and Blast2GO (Conesa et al., 2005) was used for the functional annotation of these transcript sequences. Sequences were searched (BLASTx; word size = 3, HSP length cutoff = 33) against the National Center for Biotechnology Information nonredundant database (subset Viridiplantae, taxa: 33,090) with e-value ≤ 1.0E−10, and the top-20 BLAST hits were reported. BLAST hits for each sequence were then mapped with the GO terms. Next, annotation was performed with the default parameter settings except the e-value threshold (e-value ≤ 1.0E−10). This process also generated the Enzyme Code numbers from the Kyoto Encyclopedia of Genes and Genomes pathway. After executing the whole functional annotation process (BLASTx, mapping, and annotation), transcript sequences were exported with sequence description and GO terms. Descriptions of transcript sequences were generated based on similarity levels (e-value and percentage of identity) to the subject genes, as determined by Blast2GO (Conesa et al., 2005).
Identification of Putative Orthologs and Their Trichome-Enriched Expression
Reciprocal BLAST was used to identify putative orthologs of S. pennellii genes in S. lycopersicum. An all-versus-all BLAST was performed between annotated proteins of S. pennellii v2.0 (Bolger et al., 2014a) and S. lycopersicum ITAG2.3 (Tomato Genome Consortium, 2012; Fernandez-Pozo et al., 2015) with the following parameters: minimum percentage identity = 70, minimum percentage query coverage = 50. Pairs with the lowest e-value for each BLASTp comparison were selected as putative orthologs. Gene expression profiles of S. pennellii putative orthologs in S. lycopersicum were obtained from Ning et al. (2015).
Gene Coexpression Network Analysis
Simple correlation analysis among expression profiles (fragments per kilobase of transcript per million mapped reads [FPKM] values in our 29 samples) of selected genes was performed using the “cor” function in the R programming language (R Core Team, 2014). Gene coexpression network analysis and identification of modules (clusters of strongly correlated genes) were performed using the R package WGCNA (Langfelder and Horvath, 2008). Normalized gene expression data (FPKM values) in our 38 samples were used as the input, and the 19,378 genes which passed our filtering criteria for minimum expression levels during “LOW” versus “HIGH” differential testing (Sopen12g004230 and Sopen12g004240 were merged) were selected for WGCNA analysis. A matrix of pairwise SRCCs between these 19,378 genes was created and then transformed into an adjacency (connection strength) matrix using “signed” network type and a soft threshold (β) value of 12. The value of β was determined based on the scale-free topology fit index (R2 > 0.85) and low mean connectivity. The weighted adjacency matrix was then converted to a topological overlap matrix (TOM), and genes were hierarchically clustered based on dissimilarity (dissTOM = 1 - TOM) of the topological overlap. Coexpression modules were defined as the branches of the clustering tree (dendrogram) using the following parameters in “cutreeDynamic” function: distM = disTOM, deepSplit = 2, minClusterSize = 50. A module eigengene (ME; representative of gene expression profile in each module) was defined as the first principle component in each module. A set of 41 modules was obtained. For each gene, a “module membership” was calculated as the correlation (SRCC) between expression profile of that gene and ME. Associations between external traits and modules were measured as determined by SRCC between trait data and MEs; “gene significance” values for each gene in a module of interest were calculated as the correlation (SRCC) between expression profile of a gene and trait data. The top-50 intramodular connections for selected genes were visualized with the program VisANT (Hu et al., 2005).
Inhibitor Study (S. pennellii Accession LA0716)
Compound leaves of S. pennellii LA0716 bearing five leaflets were used for the inhibitor study. Imazapyr (Sigma-Aldrich) was dissolved in distilled water, and administered as 0.1 mM and 1 mM solutions for 36 h. Cuttings were transferred to small trays containing imazapyr or control solution (deionized water) and placed under growth chamber conditions (24°C day/20°C night temperature, 150 mMol m−2 s−1 photosynthetically active radiation, 16-h photoperiod, and 75% humidity). Total RNA was extracted as described before, and RNA-Seq libraries were prepared using the TruSeq Stranded RNA LT Kit (Illumina). Libraries were sequenced on the HiSeq 4000 (Illumina) 150-bp × 150-bp paired-end sequencing platform at the Texas A & M Genomics and Bioinformatics Service Center, College Station. Raw reads were processed with the software Trimmomatic (Bolger et al., 2014b), and read-mapping to the S. pennellii genome was performed with the software Tophat2 (Kim et al., 2013) using the same parameters as described in the Quality Control and Mapping of Sequencing of Reads section, except with the following change: -library-type = fr-firststrand. Uniquely-mapped reads from TopHat2 were counted with HTSeq-Count (Anders et al., 2015) using the following options: -f bam, -r name, -s reverse, -m union, -a 20. DEGs (FDR < 0.05; log2FC > 1) were identified using the software edgeR (Robinson et al., 2010) with tagwise dispersions (prior.df = 12), and genes with >1 CPM in at least three samples were retained for the analysis. Summary results of the inhibitor study are given in Supplemental Data Set 6.
RT Quantitative PCR (S. pennellii Accession LA0716)
Stems of 10-week–old S. pennellii LA0716 plants were cut into small pieces and immediately frozen with liquid nitrogen. Stem trichomes were isolated by gentle scraping with a scalpel. Total RNA was extracted using the RNAqueous Total RNA Isolation kit (Thermo Fisher Scientific), followed by genomic DNA digestion using the TURBO DNA-free kit (Thermo Fisher Scientific), and 1 μg of total RNA was reverse-transcribed using the SuperScript IV VILO Master Mix (Thermo Fisher Scientific), following the manufacturer’s guidelines. Quantitative PCR (qPCR) was performed using SYBR Green Master Mix (Bio-Rad) with QuantStudio 6 Flex System, and the following PCR program was used: 50°C for 2 min, 95°C for 10 min, followed by 40 cycles of 95°C for 15 s and 60°C for 1 min. Differential gene expression was measured using the comparative CT (ΔΔCT) method. Three housekeeping genes (actin, ubiquitin, and TATA-box binding protein) were assessed for their use as endogenous control(s), and final normalization was performed with the ubiquitin gene Sopen02g027600. For VIGS RT-qPCR analysis, similar-sized young leaves were used for RNA extraction, and RT-qPCR primers were designed outside the VIGS-targeted region. A list of primers used in this study is given in Supplemental Table 6.
VIGS and Metabolite Profiling (S. pennellii Accession LA0716)
VIGS constructs were designed using the SGN online tool (http://vigs.solgenomics.net/). Approximately 495-bp to 510-bp target fragments were cloned into the pTRV2-LIC vector (Dong et al., 2007). Agrobacterium tumefaciens strain GV3101 carrying pTRV1 and pTRV2 vectors were grown overnight at 28°C. Overnight cultures were centrifuged at 8,000g for 5 min at 4°C, and cells were resuspended in infiltration medium (10 mM of MES at pH 5.5, 10 mM of MgCl2, and 200 μM of acetosyringone). TRV1 and TRV2 suspensions were incubated at room temperature for 3 h, and then mixed in a 1:1 ratio with a final OD600 of 1. Infiltration of S. pennellii LA0716 was performed at the first true leaf stage, and all plants were grown for about six weeks in the same growth chamber mentioned before. A phytoene desaturase gene (Sopen03g041530) was used as a positive control to assess the efficiency of VIGS (Supplemental Figure 12C).
For the analysis of acylsugar acyl chain composition, secreted acylsugars were extracted from similar-sized young leaves by dipping them in ethanol. Ethanol was evaporated, and acylsugars were dissolved in n-heptane. Transesterification reaction was performed by adding 500 μL of 20% sodium ethoxide to 1 mL of heptane mixture at room temperature for 10 min with gentle vortexing. The reaction mixture was washed twice with 500 μL of saturated sodium chloride solution, and the heptane layer was used for gas chromatography-mass spectrometry (GC–MS) analysis using Trace GC Ultra with DSQII system (Thermo Fisher Scientific). A 30-m × 0.25-mm column with 0.25-μm thickness (DB-5MS; Agilent) was used as the GC column. A quantity of 1 μL of heptane extract was used for split injection (1:50) with injector temperature of 225°C, and the following GC temperature program was used: starting temperature of 30°C held for 2 min, increased at 15°C/min up to 300°C and held for 5 min. The carrier gas flow was 1.5 mL/min. Thermo Xcalibur (v3.0.63) and Qual Browser applications were used for data acquisition and analysis. Summary of statistical tests is given in Supplemental Data Set 7.
Phylogenetic Analyses
Predicted transit peptide sequences were removed from protein sequences before multiple sequence alignment using MUSCLE (Edgar, 2004). The evolutionary history was inferred by using the maximum likelihood method based on the JTT+G model (Jones et al., 1992; lowest Bayesian Information Criterion value). Bootstrap values were obtained from 1,000 replicates. All positions containing gaps and missing data were eliminated. Evolutionary analyses were conducted in the software MEGA7 (Kumar et al., 2016). Trees were drawn to scale, with branch lengths measured in the number of substitutions per site.
Accession Numbers
RNA-Seq reads used in this study were submitted to the National Center for Biotechnology Information Sequence Read Archive under accession numbers SRP136022 and SRP136034 (imazapyr study).
Supplemental Data
Supplemental Figure 1. Previous model of acylglucose biosynthesis.
Supplemental Figure 2. Enriched GO terms associated with the 1,087 common DEGs.
Supplemental Figure 3. Differential gene expression analyses between low- and high-acylsugar–producing accessions.
Supplemental Figure 4. Expression levels of 10 selected “housekeeping genes” in different biological groups.
Supplemental Figure 5. Biosynthesis of branched-chain and straight-chain acyl molecules.
Supplemental Figure 6. Single transcript produced from the locus containing Sopen12g004230-Sopen12g004240.
Supplemental Figure 7. Correlation among expression profiles of putative defense genes.
Supplemental Figure 8. WGCNA (Langfelder and Horvath, 2008).
Supplemental Figure 9. Selected o-methyltransferase genes on chromosome 6 involved in flavonoid metabolism in S. pennellii and S. lycopersicum.
Supplemental Figure 10. Enrichment of GO terms associated with WGCNA “green” module genes.
Supplemental Figure 11. Molecular phylogenetic analysis of KAS sequences from different solanaceous species.
Supplemental Figure 12. VIGS of two candidate genes.
Supplemental Table 1. Amount of acylsugars produced by different accessions of S. pennellii.
Supplemental Table 2. Read mapping (TopHat2 alignment) results with different accessions of S. pennellii.
Supplemental Table 3. Sequence similarity between Cuphea KAS IV/KAS II-like enzymes and Sopen12g004230-Sopen12g004240.
Supplemental Table 4. Trichome-enriched expression of candidate genes’ putative ortholog in S. lycopersicum.
Supplemental Table 5. Effect of imazapyr treatment on expression levels of ABC transporter genes.
Supplemental Table 6. List of primers used in this study.
Supplemental Data Set 1. Expression levels of all S. pennellii-annotated genes and results of differential gene expression analysis between low- and high-acylsugar–producing accessions.
Supplemental Data Set 2. SRCC among selected genes’ expression profiles in 29 biological samples (different accessions of S. pennellii).
Supplemental Data Set 3. Reciprocal best hits (RBH) between S. pennellii and S. lycopersicum annotated protein sequences.
Supplemental Data Set 4. Expression of selected genes in response to imazapyr treatment.
Supplemental Data Set 5. Results of WGCNA.
Supplemental Data Set 6. Summary of the imazapyr treatment study.
Supplemental Data Set 7. Summary of statistical tests in the VIGS study.
Supplemental File 1. Alignment of KAS III sequences from different solanaceous plants.
Supplemental File 2. Alignment of KAS IV/II-like sequences from different solanaceous plants.
Supplemental File 3. Phylogenetic relationships among KAS III sequences from different solanaceous plants (mtsx file).
Supplemental File 4. Phylogenetic relationships among KAS IV/II-like sequences from different solanaceous plants (mtsx file).
DIVE Curated Terms
The following phenotypic, genotypic, and functional terms are of significance to the work described in this paper:
Acknowledgments
We thank Charlie Johnson and the staff of the Texas A&M Genomics and Bioinformatics center for performing Illumina sequencing, Rodolfo Aramayo and Ricardo Perez for technical support on the BioGalaxy server, and Yohannes Rezenom from the Texas A&M Chemistry Mass Spectroscopy Center for his assistance with GC–MS analysis. Most of the sequence analysis was performed at the Texas A&M University High Performance Research Computer, where Michael Dickens provided invaluable assistance. We sincerely appreciate Alan Pepper’s (Department of Biology, Texas A&M University) assistance with multiple aspects of this project, from initial design to critical review of the article. This work was partially supported by the US Department of Agriculture (grant 2011-38821-30891).
AUTHOR CONTRIBUTIONS
S.M. designed experiments, prepared samples, analyzed data, and wrote the article; W.J. designed experiments, prepared samples, and reviewed the article; T.D.M. designed experiments, analyzed data, and wrote the article; all authors read and approved the final article.
Footnotes
↵[OPEN] Articles can be viewed without a subscription.
- Received July 19, 2019.
- Revised September 30, 2019.
- Accepted October 16, 2019.
- Published October 18, 2019.
References
