Tomato (Solanum lycopersicum) cultivar Ailsa Craig (AC) was used in this study (Horticulture Research International, Warwick, United Kingdom). The seeds were sterilized with 10% NaClO (v/v) for 15 min, washed thoroughly with sterile water, and soaking in sterilized water overnight. After that, the seeds were sown on Petri dishes containing 1/5 strength Hoagland nutrient solution (pH 5.5). The nutrient solution consisted of KNO₃ (1.0 mM), Ca(NO₃)₂ (1.0 mM), MgSO₄ (0.4 mM), and (NH₄)H₂PO₄ (0.2 mM) and the micronutrients NaFeEDTA (20 μM), H₃BO₃ (3.0 μM), MnCl₂ (0.5 μM), CuSO₄ (0.2 μM), ZnSO₄ (0.4 μM), and (NH₄)₆Mo₇O₂₄ (1 μM), with 0.8% agar (Sigma-Aldrich, Saint Louis, United States). Perti dishes were kept in the dark at 4°C for 2 days and then germinated in a plant growth room with a daytime 16 h/24°C and 8 h/22°C night regime. Germinated seedlings with uniform primary root length (4 cm) were transferred to 1/5 Hoagland nutrient solution (pH 5.5) with (NH₄)H₂PO₄ concentration of 10 μM.

The Hidden Markov Model (HMM) file corresponding to the AMP-bind domain (PF00501) was downloaded from the Pfam protein family database¹ (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² (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³ (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.⁴

The sequences of 53 identified tomato AAEs and Arabidopsis all 60 AAE sequences according to two studies previously reported (Shockey et al., 2002, 2003; De Azevedo Souza et al., 2008) were used to create multiple protein sequence alignments using ClustalW in MEGA 7.0⁵ (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 (Yu et al., 2017).

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⁶ (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).

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.

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₄H₂PO₄) containing 5 μM Al or 5 μM of CdCl₂, or LaCl₃, or 3 μM of CuCl₂ 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.

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⁷ to identify their related cis-elements (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).

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₄H₂PO₄) 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).

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.

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).

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 data (Shockey et al., 2002, 2003; De Azevedo Souza et al., 2008) by using neighbor-joining method. In consistent with previous result (Shockey et al., 2003), AAE superfamily could be separated into six distinct subfamilies (Figure 1). Among the 53 SlAAE proteins, 15 belong to clade I, 4 to clade II, 19 to clade IV (the largest clade), 14 to clade V, and 1 to clade VI (the smallest clade). Notably, clade III, a special subfamily, contains only 19 AtAAEs (Figure 1).

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 eukaryote-type 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 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).

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 exon-intron 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.

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 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.

According to a previous study (Holub, 2001), 200 kb of a chromosomal area containing two or more genes is defined as a tandem duplication event. As shown in Figure 5, four gene pairs present as tandem duplication (T) observed on three chromosomes, i.e., T1 (Solyc02g081350, Solyc02g081360, and Solyc02g081370) and T2 (Solyc02g082870 and Solyc02g082880) on Chr02, T3 (Solyc07g043630, Solyc07g043640, Solyc07g043650, and Solyc07g043660) on Chr07, and T4 (Solyc08g075800 and Solyc08g075810) on Chr08. Therefore, in-tandem AAE duplicates comprise 21% of the whole tomato AAE superfamily.

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).

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 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² = 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.

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.

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.