Genome-Wide Identification and Gene Expression Analysis of Acyl-Activating Enzymes Superfamily in Tomato (Solanum lycopersicum) Under Aluminum Stress

In response to changing environments, plants regulate gene expression and subsequent metabolism to acclimate and survive. A superfamily of acyl-activating enzymes (AAEs) has been observed in every class of creatures on planet. Some of plant AAE genes have been identified and functionally characterized to be involved in growth, development, biotic, and abiotic stresses via mediating diverse metabolic pathways. However, less information is available about AAEs superfamily in tomato (Solanum lycopersicum), the highest value fruit and vegetable crop globally. In this study, we aimed to identify tomato AAEs superfamily and investigate potential functions with respect to aluminum (Al) stress that represents one of the major factors limiting crop productivity on acid soils worldwide. Fifty-three AAE genes of tomato were identified and named on the basis of phylogenetic relationships between Arabidopsis and tomato. The phylogenetic analysis showed that AAEs could be classified into six clades; however, clade III contains no AAE genes of tomato. Synteny analyses revealed tomato vegetable paralogs and Arabidopsis orthologs. The RNA-seq and quantitative reverse-transcriptase PCR (qRT-PCR) analysis indicated that 9 out of 53 AAEs genes were significantly up- or downregulated by Al stress. Numerous cis-acting elements implicated in biotic and abiotic stresses were detected in the promoter regions of SlAAEs. As the most abundantly expressed gene in root apex and highly induced by Al, there are many potential STOP1 cis-acting elements present in the promoter of SlAAE3-1, and its expression in root apex was specific to Al. Finally, transgenic tobacco lines overexpressing SlAAE3-1 displayed increased tolerance to Al. Altogether, our results pave the way for further studies on the functional characterization of SlAAE genes in tomato with a wish of improvement in tomato crop in the future.


INTRODUCTION
Aluminum (Al) toxicity is one of the major limiting factors affecting the crop productivity in acidic soils, which occupy nearly 50% of the potential arable lands of the world (von Uexküll and Mutert, 1995). When soil pH is lower than 5.5, ionic Al, mainly Al 3+ , predominates in soil solution, which is highly toxic to plants. The initial and most visible symptom of Al toxicity is inhibition of root elongation by ravaging cell structure of the root apex and thus limiting the mineral nutrient and water uptake and, consequently, hindering the plant growth and development (Kochian, 1995;Ma, 2007;Ryan et al., 2011;Liu et al., 2014). To adapt to Al-toxic environment, plants have evolved two major types of Al-tolerance mechanisms, namely, external exclusion (preventing Al from entering cells of root apex) and internal tolerance mechanisms (detoxifying Al via complexation and sequestration) (Kochian et al., 2004(Kochian et al., , 2015Liu et al., 2014). Substantial advances have been made toward elucidating the physiological and molecular mechanisms by which plants cope with Al stress (Kochian et al., 2015;Yang et al., 2019). Recently, it has been shown that metabolic change might play important roles in response to Al stress (Lou et al., 2016a,b;Xian et al., 2020). However, the molecular basis of the role of metabolic alterations in Al stress response still needs further elucidation.
The activation of carboxylic acids provides the precursors for pathways that lead to the metabolism of a diverse variety of metabolites, including lipids, amino acids (aa), sugars, and secondary metabolites. In plants, the acyl-activating enzymes (AAEs) superfamily consists of acyl-coenzyme A synthetases (ACSs), 4-coumarate:coenzyme A ligases (4CLs), luciferases, and non-ribosomal peptide synthetases, which are involved in many primary and secondary metabolic pathways. All members of the AAE family have low sequence similarity to each other but share many highly conserved motifs, such as the AMP-binding domain (Shockey and Browse, 2011). As the members of the Arabidopsis AAE superfamily were systematically analyzed and identified, more and more metabolic functions of plant AAE superfamily genes have been reported. For example, the Arabidopsis peroxisomal-localized OPCL1 (OPDA-CoA ligase) gene, At1g20510, is involved in the biosynthesis of jasmonic acid in Arabidopsis (Koo et al., 2006;Kienow et al., 2008). The rice fatty acyl-CoA synthase gene, OsACOS12, is involved in regulating lipid metabolismmediated tapetum-programmed cell death, which ultimately affects the male fertility of rice (Yang et al., 2017). The petunia malonyl-CoA synthase gene, PhAAE13, is specifically involved in anthocyanin biosynthesis in flowers . Recently, a rice 4CL4 belonging to 4-coumarate:coenzyme A ligases was reported to be involved in Al resistance (Liu et al., 2020).
Al-induced secretion of organic acids, including citrate, malate, and oxalate, has been well-documented as a very important mechanism by which plants resist the Al toxicity (Yang et al., 2019). Although transporters responsible for Alinduced citrate and malate secretion, respectively, have been characterized in a variety of plant species, genes encoding oxalate transporter remain unclear (Yang et al., 2019). Accumulating evidence suggests that oxalic acid has an important role in plant responses to both biotic (Molano-Flores, 2001;Jang et al., 2016) and abiotic stresses, including calcium regulation, ion homeostasis, metal stress, and other pathways (Palmieri et al., 2019). It has been reported that one of the AAE family members, AAE3 (acyl-activating enzyme3), is involved in oxalic acid metabolism-mediated plant growth and development and in resistance to biotic and abiotic stresses. For example, Foster et al. (2012) identified an AAE3 gene encoding an oxalyl-CoA synthase in Arabidopsis and found that it participated in seed development and fungal pathogen defense by catalyzing CoAdependent oxalate metabolism (Foster et al., 2012). Lou et al. (2016a) found that rice bean (Vigna umbellata) VuAAE3 is involved in oxalate degradation and Al tolerance. Recently, Xian et al. (2020) also found that wild soybean GsAAE3 similarly influences its tolerance to Cd and Al stress by catalyzing the oxalate metabolism. Therefore, it appears that AAE family proteins might have important roles in Al stress responses by regulating metabolic pathways.
In this study, we identified 53 AAE genes from tomato genome and found 9 differentially expressed AAE genes, including SlAAE3-1 and SlAAE3-2, under Al stress. The expression pattern analysis of SlAAE3-1 suggested that its expression is specific to Al stress. Therefore, our results contribute not only to enrich the molecular mechanism of Al stress response in tomato but also to provide a theoretical basis for improving tomato Al tolerance through genetic improvement and molecular breeding techniques.

Identification of Acyl-Activating Enzyme Superfamily in Tomato
The Hidden Markov Model (HMM) file corresponding to the AMP-bind domain (PF00501) was downloaded from the Pfam protein family database 1 (El-Gebali et al., 2019). HMMER 3.2 was used to search against the AAE superfamily genes from the annotated tomato genome obtained from Phytozome version 12.1 2 (Finn et al., 2011;Goodstein et al., 2012). All candidate genes that may contain AMP-binding domain based on HMMER results were further examined by confirming the existence of the AMP-binding core sequence using the PFAM and the SMART program 3 (Letunic and Bork, 2018). The length of aa sequences, protein molecular weights (MWs), and isoelectric point of identified tomato AAE superfamily proteins were obtained by using tools from the ExPasy website. 4

Phylogenetic Analysis of Acyl-Activating Enzyme Supfamily Members
The sequences of 53 identified tomato AAEs and Arabidopsis all 60 AAE sequences according to two studies previously reported (Shockey et al., 2002(Shockey et al., , 2003De Azevedo Souza et al., 2008) were used to create multiple protein sequence alignments using ClustalW in MEGA 7.0 5 (Kumar et al., 2016) with default parameters. The alignment results were used to construct a phylogenetic tree using the neighbor-joining method with 1,000 bootstrap replicates. The phylogenetic tree was displayed using the R package ggtree .

Gene Structure and Conserved Motif Analysis
The exon-intron distribution of each tomato AAE superfamily genes (SlAAEs) was analyzed by comparing predicted coding sequences with their corresponding genomic sequences using TBtools program (Chen et al., 2020). Conserved motifs of tomato AAE protein sequences were investigated using the online software MEME5.0.4 6 (Bailey et al., 2009) with the following motif parameters: number of repetitions (any), maximum number of motif (20), and the optimum width of each motif (between 6 and 100 residues).

Chromosomal Distribution and Gene Duplication Analysis
All SlAAEs were mapped to 12 tomato chromosomes based on physical location information from the database of tomato genome using TBtools (Chen et al., 2020). Multiple Collinearity Scan Toolkit, McScanX (Wang et al., 2012) with the default parameters was used to analyze the tandem repeats and segmental duplication events of SlAAEs superfamily in the tomato genome and synteny of AAE superfamily genes between tomato and Arabidopsis.

Expression Analysis of Aluminum-Responsive SlAAE Genes
To investigate Al-responsive SlAAEs, the analysis of RNA-seq data (Jin et al., 2020) and qRT-PCR were performed. For qRT-PCR, seedlings were subjected to the modified 1/5 Hoagland nutrient solution (pH 5.0; 10 µM NH 4 H 2 PO 4 ) containing 5 µM Al or 5 µM of CdCl 2 , or LaCl 3 , or 3 µM of CuCl 2 for 6 h. RNA samples were extracted from both root tips (1 cm in length) after treatment. One microgram of DNA-free RNA was transcribed into first strand cDNA by PrimeScript RT Master Mix (TaKaRa). The qRT-PCR was carried out with the Roche LightCyler 480 instrument using SYBR Green chemistry (Toyobo, Osaka, Japan). The reaction conditions were 40 cycles at 95 • C for 15 s, 60 • C for 10 s, and 72 • C for 15 s. The primer sequences used in this study are listed in Supplementary Table 1. Expression data of target genes were normalized with tomato GAPDH (Wang et al., 2016) by the Ct method. Each reaction was performed with three repeats from different biological samples.

Promoter Analysis
The promoter data were obtained from Phytozome version 12.1 (see text footnote 2). The promoter analysis was conducted by searching 2.0 kb upstream sequences of the coding sequences against the PlantCARE database 7 to identify their related ciselements (Lescot et al., 2002). After sorting the cis-elements obtained from PlantCARE, the results were visualized and mapped to the AAE promoter using the TBtools software (Chen et al., 2020).

Overexpression of SlAAE3-1 in Tobacco and Aluminum Tolerance Evaluation
The open read frame of SlAAE3-1 was amplified by PCR using gene-specific primer pair TGGTACCATGGAGAGTATGACGCTC and CCGGGATCCCTACGCTCCAAATTTAGG, cloned into pCAMBIA1300 vector driven by Cauliflower mosaic virus 35S promoter and transformed into Agrobacterium tumefaciens (strain GV1301). Tobacco plants were transformed as described by Horsch et al. (1985). Transgenic lines carrying SlAAE3-1 were selected using PCR with the primers described above. For evaluating Al tolerance of SlAAE3-1-overexpressing lines, seeds from T2 homozygous and wild-type lines were first sterilized, soaked, and germinated as described above. When the length of the primary root had reached about 4 cm, the seedlings were transferred to the modified 1/5 Hoagland nutrient solution (pH 5.0; 10 µM NH 4 H 2 PO 4 ) containing 5 µM Al for 4 days. The solution was renewed every 2 days. The Al sensitivity was evaluated by relative root elongation expressed as (root elongation with Al treatment/root elongation without Al) × 100. After the treatment, root apex was stained with propidium iodide (PI) solution (5 µg/ml) for 1 min, washed with deionized water for 30 s, then observed, and captured using confocal laser scanning microscopy (Zeiss LSM710, Jena, Germany).

Statistical Analysis
The Student's t-test was performed in Microsoft Excel (version 2016, Microsoft Corp., Redmond, WA, United States). Data are given as means ± standard deviation (SD) of three independent biological replicates. A p-value < 0.05 was considered to be statistically significant.

Identification of the Acyl-Activating Enzymes Superfamily in Tomato
We identified 53 members of AAEs superfamily in S. lycopersicum by searching the AAE consensus motif (PF00501) equal to PROSITE PS00455 (Shockey et al., 2003) following a previously described analysis pipeline (Jin et al., 2020). Gene characteristics, including the length of the coding sequence, the length of the protein sequence, the protein MW, isoelectric point (pI), and protein sequence, were analyzed (Supplementary Table 2). Among the 53 SlAAE proteins, Solyc07g043650 was identified to be the smallest protein with 455 aa, whereas the largest one was Solyc01g006640 (2,320 aa). The MW of SlAAE proteins ranged from 50.4 to 256.8 kDa, and the pI varied from 5.13 (Solyc08g076300) to 9.43 (Solyc03g005090).

Phylogenetic Analysis and Classification of SlAAE Genes
To probe the phylogenetic relationships among these 53 SlAAEs, we constructed an unrooted phylogenetic tree for SlAAEs together with 60 AtAAEs retrieved from previously published FIGURE 1 | The phylogenetic analysis of tomato (S. lycopersicum) acyl-activating enzymes (AAEs) (SlAAEs). The phylogenetic analysis of AAEs from tomato and Arabidopsis using the complete protein sequences. The neighbor-joining (NJ) tree was constructed using the MEGA 7.0 software with the pairwise deletion option, and 1,000 bootstrap replicates were used to assess tree reliability. AAEs from tomato and Arabidopsis fell in six separate subfamilies as I-VI.
According to the feature of known long-chain acyl-CoA synthetase (LACS) proteins (Iijima et al., 1996;Fulda et al., 1997), 11 SlAAE members from clade I contain the eukaryotetype linker domain, a motif of between 30 and 70 aa residues (Figure 2), and may be active against long-chain fatty acids. However, similar to previously characterized AtLACSs (Shockey et al., 2002), the other four members (i.e., SlAAE3-1, SlAAE3-2, Solyc02g069920, and Solyc09g092450) probably do not produce the activity of the LACS enzyme even though they showed highly sequence similarity to 11 members. So far, the biological functions of 13 AAEs from clade I remain unknown except Solyc01g079240 (SlLACS1) and Solyc01g109180 (SlLACS2), both of which are reported to be involved in wound-induced suberization of tomato fruit (Han et al., 2018).
Six members of clade II could be further separated into two subgroups. Solyc02g094640 and Solyc07g017860 from clade II share 78% sequence similarity with AT5G63880 that was reported functioning as acetyl-CoA synthetase involved in lipid synthesis in seeds (Ke et al., 2000). The results of multiple sequence alignment among Solyc03g114460, Solyc12g038400, AAE17, and AAE18 revealed 69% sequence identity. However, all proteins from clade II share only 51% sequence similarity on average. These results suggest that these two subgroups might have distinct functions FIGURE 4 | The sequence information of each conserved motif (1-15). Conserved motifs of acyl-activating enzyme (AAE) superfamily members were found using the MEME online tool with parameters: number of repetitions (any), maximum number of motif (20), and the optimum width of each motif (between 6 and 100 residues).
Frontiers in Plant Science | www.frontiersin.org and require further biochemical assays to elucidate their functions.
Clade III consisted of 19 Arabidopsis AAE family members, termed adenylases, which are considered to participate in multiple important plant hormones (e.g., JA, IAA, and SA) signaling pathways through ATP-dependent adenylation of these hormones (Staswick et al., 2002;Shockey et al., 2003). Notably, no tomato AAEs were included in this clade. Among 19 tomato AAE genes in clade IV, only Solyc07g043630 (SlAACS1) has been reported to be involved in biosynthesis of acylsugars in tomato trichomes (Fan et al., 2020), while the function of the remaining genes has not been characterized yet.
Clade V contains 13 putative Arabidopsis 4-coumarate CoA ligases (At4Cls) (Shockey et al., 2003) and 14 putative tomato 4CLs. The 4CLs play a vital role in enhancing the mechanical support of plants and protecting plants from biotic and abiotic stresses depended on biosynthesis of lignins, flavonoids, and other compounds (Lavhale et al., 2018). For example, 4CL gene involving in biosynthesis of lignin is upregulated in tomato (S. lycopersicum) upon Alternaria solani inoculation (Shinde et al., 2017). In rice, Al repressed the expression of 4CL4, resulting in less lignin accumulation and more 4-coumaric acid and ferulic acid accumulation (Liu et al., 2020). In Arabidopsis, At4CL1/2 required for biosynthesis of lignin was upregulated, while At4CL3 involved in flavonoid was downregulated by wounding (Lee and Douglas, 1996;Ehlting et al., 1999;Soltani et al., 2006).

Gene Structure and Protein Motif Analysis of SlAAE Genes
During the evolution of multigene families, diversification of gene structure may facilitate evolutionary co-option of genes for new functions to adapt to changing environments (Shang et al., 2013;Hu et al., 2015). To better understand the structural diversity of SlAAE genes, the conserved motifs and exonintron organizations were analyzed according to the phylogenetic relationships (Figure 3).
The potential conserved motifs of all SlAAE proteins were presented according to the MEME motif analysis as described (Xie et al., 2018). As a result, 15 different motifs were identified from all SlAAEs, which were successively named as motifs 1-15 ( Figure 3A). The sequence information of each motif is displayed in Figure 4. Among these motifs, motifs 3, 5, and 9 belong to the AMP-binding domain, which are widely distributed on all SlAAEs. As expected, SlAAEs members within the same cluster in the phylogenetic tree commonly share a similar motif compositions (Figure 3A), indicating that they might have functional similarities.
As shown in Figure 3B and Supplementary Table 2, SlAAEs genes contain 0-22 introns. For instance, while Solyc02g037490 and Solyc04g080480 have no intron, Solyc01g099100, Solyc09g075770, and Solyc09g075790 have 22 introns. Others have at least one intron, but genes with 8, 12, 14, 15, 20, and 21 introns were not observed. Most of SlAAE members in the same group have a similar gene structure. For example, most of group IV members contained only a single intron. These results suggest that SlAAEs possessing similar gene structures and motifs were clustered in the same group and might have evolved similar functions in tomato.

Chromosomal Localization of SlAAE Genes
To probe the chromosomal distribution of all SlAAE genes, all members were mapped to tomato chromosomes based on physical location information derived from the database of FIGURE 6 | Syntenic analysis of AAE genes in the tomato genome. (A) Identification of paralog pairs in syntenic blocks within the tomato genome. Ten paralog pairs were identified. (B) Synteny analysis of AAE genes between tomato and Arabidopsis. Three syntenic pairs were found in this analysis. Gray lines in the background indicate the collinear blocks within tomato and Arabidopsis, while the red lines highlight the syntenic AAE gene pairs. tomato genome. The result showed that all of the 53 SlAAEs could be mapped onto 10 out of 12 tomato chromosomes (except Chr05 and Chr10) in an increasing order from short-arm to long-arm telomere (Figure 5), and most of SlAAE genes were distributed in the chromosome ends. In addition, bias change in gene number was inspected. Among 10 chromosomes, Chr02 contained largest number of AAE genes, while Chr04 had least.

Syntenic Analysis of AAE Genes in Tomato Genome
Besides the tandem duplication events, we also investigated the segmental duplication in the tomato genome relating to the recurring polyploidization events, which generated gene duplicates that have usually been retained in extant tomato genome (Wang and Paterson, 2011). In this study, Eight pairs (named pair 1-8) of syntenic AAE paralogs were observed within the tomato genome ( Figure 6A). According to the phylogenetic tree (Figure 1), the tomato paralogs belong to clade I (syntenic pairs 1 and 8), clade IV (syntenic pairs 2, 3, and 4), as well as clade V (syntenic pairs 5, 6, and 7), thus allowing us to propose the conserved functions between the syntenic pairs. Furthermore, we constructed a comparative syntenic map between tomato and Arabidopsis ( Figure 6B). Three syntenic pairs comprising four SlAAE genes were identified in syntenic blocks. We also found that duplicated AAE genes exhibited conserved synteny with Arabidopsis genes. For example, Solyc02g082870 is microsyntenic to Solyc03g032210 (paralog pair 4), both of which are syntenic to AT2G17650 (ortholog pair 2, Figure 6B). Solyc08g082280 is microsyntenic to Solyc01g095750 (paralog pair 1) and syntenic to three Arabidopsis genes (ortholog pair 3). Based on these information, we could have inferred the functions of these tomato genes based on the functions of Arabidopsis genes, though the functions of these genes in Arabidopsis are yet to be investigated. Nevertheless, the biological functions of these ortholog genes could have been conservatively evolved since the last common ancestor of tomato and Arabidopsis, which is estimated to have existed approximately 150 million years ago (Ku et al., 2000).

Expression Profiles of SlAAEs Under Aluminum Stress
We have previously demonstrated that rice bean (V. umbellata) VuAAE3, a gene showing a high sequence identity to SlAAE3-1 (Solyc03g025720) and SlAAE3-2 (Solyc06g035960), was involved Expression specificity of SlAAE3-1 gene. Seedlings were subjected to 1/5 Hoagland solution containing 5 µM AlCl 3 , or 5 µM CdCl 2 , 5 µM LaCl 3 , or 3 µM CuCl 2 for 6 h. Data are the means ± SD (n = 3); the asterisk indicates significant differences between control and treatment at P < 0.001 using one-way ANOVA. (C) Correlation of gene expression levels between the RNA-seq data and quantitative reverse-transcriptase (qRT-PCR) analysis. Fourteen SlAAEs potentially responding to Al were selected and subjected to the qRT-PCR analysis using the same RNA as for RNA-seq. Both x-and y-axes are shown in Log 2 scale.
in Al tolerance (Lou et al., 2016a). This prompted us to investigate the potential function of SlAAE genes in Al stress response after a systemic analysis of the SlAAE gene family. On the basis of tomato root tip Al stress-responsive expressed genes identified from the results of RNA-seq (SRP227103) (Jin et al., 2020), we found that 50 out of 53 SlAAE genes could be detected by RNA-seq in the tomato root apex (Supplementary Table 3). However, six AAE genes showed FPKM value lower than six and were hardly expressed either without or with Al stress. In addition, the expression of 35 AAE genes was not induced by Al stress in tomato root apexes (Supplementary Table 3). Finally, only 9 out of 53 SlAAEs were identified to be differentially regulated, 8 were upregulated, and 1 were downregulated by 5 µM of Al (Supplementary Table 4 and Figure 7A). Notably, Solyc03g025720 (SlAAE3-1) was highly expressed and greatly induced by Al compared with others. We then examined the specificity of SlAAE3-1 expression by exposing tomato seedlings to various metals, including Al, Cd, La, and Cu. The expression of SlAAE3-1 was greatly induced by Al but not by other metals (Figure 7B). To verify the reliability of the RNA-seq data, 14 SlAAE genes were selected for the qRT-PCR analysis. As shown in Figure 7C, all 14 SlAAE genes displayed similar expression patterns to that obtained using RNA-seq. A good correlation (R 2 = 0.8415) was observed for their expression in plot qRT-PCR results against that of RNA-seq, indicating that the RNA-seq data accurately reflected the transcriptional changes induced by Al stress.
FIGURE 8 | Phylogenetic and cis-element analysis of SlAAE family promoters. A total of 53 of promoter sequences from tomato genome were scanned using PlantCARE. Sorted cis-elements were then mapped on promoters of the corresponding AAEs and visualized using TBtools. Colored rectangles represent different cis-element.

Identification of Stress-Responsive cis-Acting Elements of SlAAEs
To explore the cis-acting elements probably implicated in the expression regulation of SlAAEs, 2.0 kb sequences upstream of the start codon of SlAAEs were analyzed using the PlantCARE database (Lescot et al., 2002). The results showed that most of the predicted cis-acting elements were associated with phytohormone responses. In addition, there were MYB-binding site, light-responsive element, 60K protein site, ATBP-1 binding site, and endosperm-specific negative expression. However, we did not find Al-responsive element predicted from the PlantCARE (Figure 8). Sensitive to proton rhizotoxicity 1 (STOP1) is a C2H2-type zinc finger transcription factor that regulates expression of many downstream genes involved in Al tolerance by binding to STOP1 cis-acting elements GGN(T/g/a/C)V(C/A/g)S(C/G) present in their promoters (Tsutsui et al., 2011;Liu et al., 2016). Therefore, we analyzed the 2 kb length promoter sequence of SlAAE3-1 and identified 28 cis-acting elements that have potentials to interact with STOP1 ( Table 1), suggesting that the STOP1 regulatory module may also be present in tomato.

Effect of SlAAE3-1 Overexpression in Tobacco on the Tolerance to Aluminum
To confirm that the identified differentially expressed SlAAE genes are exactly involved in tolerance to stresses, we developed transgenic tobacco lines overexpressing SlAAE3-1 that is most significantly induced by Al in the tomato root apex. Three independent SlAAE3-1 overexpressing tobacco lines (i.e., OE1, OE2, and OE3) were selected for examining their tolerance to Al stress. Under normal condition, the wild-type (WT) and transgenic lines showed no difference in root elongation. However, in the presence of 5 µM of Al, the elongation of the primary root of transgenic lines was significantly greater than WT lines (Figures 9A,B). In addition, the PI staining was used to check cell damage. In the absence of Al, PI was hardly stained both in the WT roots and transgenic lines ( Figure 9C). However, Al stress resulted in the red fluorescence signals to be more severely accumulated in the WT root apex than in the transgenic lines (Figure 9C), suggesting that transgenic tobacco lines were more tolerant to Al stress than WT plants. Therefore, SlAAE3-1 plays roles in the tolerance to Al stress, consisting with its transcriptional regulation by Al.

DISCUSSION
In higher plants, AAEs, also called acyl adenylate-forming (Conti et al., 1996;Chang et al., 1997) or AMP-binding proteins (Fulda et al., 1997), are involved in numerous metabolic pathways, such as fatty acid β-oxidation (Schnurr et al., 2004), oxalate catabolism (Foster et al., 2016), and malonic acid degradation (Chen et al., 2011). In the current study, we systemically analyzed the AAE superfamily in tomato and identified a total of 53 members (Supplementary Table 2). We further divided tomato AAE superfamily into five distinct clades based on the phylogenetic TABLE 1 | The number and position of sensitive to proton rhizotoxicity 1 (STOP1) cis-acting elements of SlAAE3-1 promoter.

Motif sequence
Position of STOP1 cis-element analysis (Figure 1). SlAAEs and AtAAEs from clades I, II, IV, V, and VI showed that these genes were not only homologous but could be evolved from a common ancestor. However, AAEs from clade III indicated that this clade of AtAAE genes had a different ancestor with tomato (Figure 1). It is also interesting to note that clade III AAEs encompass 19 Arabidopsis plant hormone adenylases including JAR1. Substrate-dependent ATP-32 P-PPi isotope exchange experiment demonstrated that JAR1 is specifically active on JA, while some members from this clade are active on auxin (Staswick et al., 2002). Here, we found that all 53 tomato AAE superfamily members had hardly been reported functionally, except SlLACS1 (Solyc01g079240) and SlLACS2 (Solyc01g109180) (Han et al., 2018). Therefore, the systematic analysis and identification of the tomato AAE superfamily in this study facilitate further studying on the biological function of the superfamily. As one of the most important metabolisms in plants, phenylpropanoid metabolism provides precursors for more than 8,000 metabolites contributing to plant development and plantenvironment interplay (Lavhale et al., 2018;Dong and Lin, 2021). The reaction catalyzed by 4-coumarate-CoA ligase (4CL) is the third step of the first three shared common steps of the general phenylpropanoid pathway, which is responsible for channelizing precursors for various phenylpropanoids (Fraser and Chapple, 2011;Lavhale et al., 2018;Dong and Lin, 2021). In plants, 4CL enzymes belong to AAE superfamily and catalyze the reaction that converts methoxy or hydroxycinnamic acid derivatives to corresponding CoA thioesters (Shockey et al., 2003;Lavhale et al., 2018). In addition, the 4CL enzymes play vital roles in plant physiology or in responses to biotic and abiotic stresses (Fritzemeier et al., 1987;Lee and Douglas, 1996;Sun et al., 2013;Abdollahi Mandoulakani et al., 2017;Blanco-Ulate et al., 2017). In rice, 4CL4-knockout mutants increase FIGURE 9 | Overexpression of SlAAE3-1 in tobacco improved Al tolerance. (A) Al tolerance phenotype of SlAAE3-1 overexpressing tobacco lines. (B) Primary root elongation. Data are the means ± standard deviation (SD) (n = 10); the asterisk indicates significant differences between control and treatment at P < 0.001 using one-way ANOVA. (C) Propidium iodide staining, fluorescence signals were analyzed using confocal microscopy, bars = 100 µm.
the Al tolerance by reducing the binding of Al to the cell walls caused by increased accumulation of 4-coumaric acid and ferulic acid that strengthens the cross-linking of the hemicellulose (Liu et al., 2020). In this study, we also identified 14 putative tomato 4CL and 4CL-like enzymes in clade V, but their biological functions have to be characterized in future. According to the RNA-seq data, we found that the expression level of 2 4CL genes, Solyc03g097030 and Solyc06g035960, increased under Al stress for 6 h ( Figure 7A). 4CLs in dicots, such as tomato and Arabidopsis, could be grouped into two clusters, namely, type I and type II (Supplementary Figure 1). Type I is mainly involved in lignin biosynthesis, whereas type II cluster is involved in phenylpropanoid biosynthesis other than lignin (Gui et al., 2011;Lavhale et al., 2018). However, in monocot plants, such as rice, five Os4CLs were mainly categorized into type III except Os4CL2 that belongs to type II (Gui et al., 2011;Sun et al., 2013). As shown in Supplementary Figure 1, Solyc03g097030 belongs to type I, indicating that Solyc03g097030 responds to Al stress possibly by regulating lignin biosynthesis. However, it remains possible that Solyc03g097030 might be involved in other biosynthesis pathways, not just lignin biosynthesis. This proposition is supported by a recent report that although Os4CL4 belongs to type III, it does participate in the regulation of lignin biosynthesis (Liu et al., 2020).
In general, the presence of multiple paralogs in multigene families may relate to the recurring polyploidization events of the angiosperm lineage, which generated gene duplicates that have often retained in extant plant genomes (Wang and Paterson, 2011). It has been shown that a genome-wide duplication event happened in tomato about 83-123 Myr (Sato et al., 2012). Over time, these gene duplicates may have culminated in sub-or neofunctionalization and, subsequently, acquired new functions that are occasionally retained, thus resulting in functional diversity and proliferation of genes derived from a common ancestor gene (Veitia, 2005). In this study, we revealed four tandem duplication segments (Figure 5), which may result in an intensification of gene expression. For example, maize with in-tandem MATE genes (three-copy allele) show a greater Al tolerance as enhanced overall expression of these genes (Maron et al., 2013). Orthologs and paralogs are two essentially different types of homologous genes that are associated with speciation or duplication (Koonin, 2005;Gabaldón and Koonin, 2013). In current study, eight pairs of syntenic AAE paralogs were found within the tomato genome (Figure 6A), and three ortholog pairs were identified in syntenic blocks between tomato and Arabidopsis ( Figure 6B). According to the phylogenetic tree, we found that tomato paralog pair 4 belong to the clade IV (Figures 1, 6A), and this paralog pair showed synteny with Arabidopsis AT2G17650 (ortholog pair 2 in Figure 6B), suggesting that these genes could share conserved functions. These results suggested that the analysis of synteny of genes contributes to inferring novel gene functions based on known genes.
We identified nine tomato SlAAE genes that rapidly responded to Al stress in the tomato root apex, among which SlAAE3-1 was most abundantly expressed and dramatically upregulated (Figure 7 and Supplementary Table 4). Previous studies have shown that rice bean VuAAE3 and wild soybean GsAAE3 are implicated in Al tolerance by regulating oxalate acetylation (Lou et al., 2016a;Xian et al., 2020). Here, we demonstrated that tomato SlAAE3-1 plays the same role with respect to Al tolerance (Figure 9). However, the role of other Alresponsive SlAAE genes in Al tolerance has to be investigated. There is considerable evidence that AAE genes play critical roles in the plant growth and development (Schnurr et al., 2004;Lü et al., 2009;Chen et al., 2011;Yang et al., 2017) and improving tolerance to abiotic stresses, including salinity , drought (Bessire et al., 2007;Zhao et al., 2019), and metal stress (e.g., Cd and Al) (Lou et al., 2016a;Xian et al., 2020). Therefore, the functional roles of these SlAAEs in various stress responses could be inferred. In accordance with this supposition, many cisacting elements related to biotic and abiotic stresses have been identified to be present in their promoters (Figure 8).
In summary, we provided the first integrated bioinformatic information of SlAAE superfamily in tomato including gene identification, structure, chromosomal location, duplication, and expression regulation by Al stress. A total of 53 SlAAE genes were identified, which is essential for the functional characterization of SlAAE3 genes in tomato in future. Furthermore, the RNA-seq data and qRT-PCR analysis have identified nine Al-responsive SlAAE genes and characterized one of them, SlAAE3-1, to be implicated in Al tolerance, which pave the way for identifying novel genes involved in Al tolerance. In addition, the central role of SlAAE members in diverse metabolisms shed light on the importance of this family in responding to both biotic and abiotic stresses.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
JY and WC conceived the research. JJ, QH, and PL performed the experiments. HL provided the technical assistance. JJ and QH analyzed the data. JJ, WC, and JY wrote the manuscript. All authors contributed to the article and approved the submitted version.

FUNDING
This work was financially supported by grants from the Zhejiang Provincial Natural Science Foundation (LY19C150006), the National Natural Science Foundation of China (31601765), the China Postdoctoral Science Foundation (to WC), and the China Scholarship Council.