The variety Micro-Tom was used as wild-type (WT) for the genetic materials in this study. Tomato plants were grown in a greenhouse at 24°C under a 16-h-light and 8-h-dark cycle at Hainan University, Haikou, China. RNA-sequencing (RNA-seq) samples were harvested at different stages during the fruit development of tomato. The skin of fruits was collected from days post-germination (DPG) at immature green (IMG in 60 DPG), mature green (MG in 80 DPG), breaker (Br in 90 DPG), and ripe fruit (Ripe in 100 DPG) stages, four periods with three biological replicates were collected. All samples were collected before the greenhouse lights were turned off at approximately 7:00–9:00 p.m. Samples were collected from the plants, immediately placed in liquid nitrogen.

Tomato ITAG3.2 genome data were obtained from Sol Genomics Network.¹ The Hidden Markov model (HMM) data of AP2 domain were downloaded from Pfam database² with domain number PF00847. We scanned the tomato genome by HMMER 3.2.1 with default parameters. The presence of AP2 sequences was confirmed by using Pfam (see text footnote 2) and Simple Modular Architecture Research Tool (SMART)³ to further check for all candidate genes that may contain AP2/ERF structural domains based on HMMER software results. We expanded ERF gene numbering based on the naming system provided by Liu et al. (2016). For genes from other subfamilies, we numbered them sequentially according to their chromosomal locations. The sequence length, molecular weight, and isoelectric point (pI) of the identified AP2/ERF family proteins were obtained using the tools on the ExPasy website.⁴ The predicted subcellular localization data for all AP2/ERF genes were obtained via the GenScript website.⁵

We collected tomato skin from the four key fruit development periods (IMG, MG, Br, and Ripe) for RNA-seq in this study. Total RNA-seq of all samples was extracted by TRIzol method. Besides, the transcriptome libraries were constructed by using an Illumina TruSeq RNA kit according to the manufacturer’s recommendations. Sequencing was performed using an Illumina HiSeq X Ten system. RNA-seq data from each sample were obtained for three biological replicates. The quality of RNA-seq data was examined using Fastp (0.20.1), and adapter sequences and low-quality reads were filtered by Fastp with default parameters. Reads were mapped to S. lycopersicum reference genome (ITAG 3.2) using Hisat2 (Kim et al., 2019). FeatureCounts (v2.0.1) was used to estimate the expression levels of all genes (Liao et al., 2019), with the command: featureCounts −T −p −t exon −g gene_id −a. Transcripts Per Million (TPM) was used to normalize the expression levels of all samples. Different expression genes (DEGs) were identified using DESeq R package (Anders and Huber, 2010), with the threshold for fold change ≥ 2 and the significance level <0.01.

Multiple sequence alignment of the AP2 structural domain of the AP2/ERF protein from Arabidopsis thaliana (Nakano et al., 2006) was performed using the software Gblocks to obtain the conserved structural sequence.⁶ The alignment was used through Molecular Evolutionary Genetics Analysis (MEGA) X using the ClustalW method with default parameters (Kumar et al., 2018). Neighbor-Joining (NJ) tree construction was performed using MEGA X with the Poisson model, pairwise deletion with 1,000 bootstrap replications.

All SlAP2/ERFs distribution information on chromosomes were based on the tomato genome database and visualized by TBtools. The results visualized using the Circos software. Tandem duplication and segment duplication were explored by using the Multiple Collinearity Scan toolkit X (MCScanX) (Wang et al., 2012). Synteny analysis was performed within the tomato genome and between multiple species also using MCScanX, and the genome data of other species were downloaded from the ensemble website.⁷

Weighted gene co-expression network analysis was performed to generate the gene expression network based on the correlation patterns among all high expression genes (Fragments Per Kilobase per Million, FPKM > = 1). The construction of a weighted co-expression network needs the soft-thresholding powers β, which were calculated by the pick soft threshold function of the R Package (WGCNA version 1.70-3) (Langfelder and Horvath, 2008). We chose the power 10 as it is the lowest power for which the scale-free topology fit index curve flattens out reaching a high value (above 0.9). The network of modules of co-expressed genes (edge weight >0.10) was represented by Cytoscape 3.7.1 (Otasek et al., 2019). The expression levels of SlAP2/ERF genes and SGA pathway genes generated in this study were used for co-expression analysis, and Pearson’s correlation coefficients (PCC) were used to measure the co-expression relationships between genes.

Whole genome sequencing (WGS) raw data were used from previous studies (Lin et al., 2014; Zhu et al., 2018). A total of 1,476,413 single nucleotide polymorphisms (SNPs) with minor allele frequency (MAF) >5% and missing rate <10% from 349 accessions were used as genotypic data. The phenotypic data used the relative content of esculeoside A from our previous study (Zhu et al., 2018), and then, we performed the genome-wide association analysis. Efficient Mixed-Model Association eXpedited (EMMAX) software was used to test the associated relationships between genotype and phenotype (Kang et al., 2010), and population structure was modeled as a random effect in Efficient Mixed-Model Association (EMMA) using Admixture software, the kinship (K) matrix was calculated by emmax-kin, the parameter is emmax-kin −v −h −s −d 10. The effective number of independent SNPs was calculated using Genetic type 1 Error Calculator (GEC) software (Li et al., 2012). The calculated genome-wide significance thresholds, based on a nominal level of 0.05, were p = 4.19 × 10^––7 for the whole population. Linkage disequilibrium (LD) was calculated based on the R² value by Plink (v1.9).⁸

The full-length coding sequence (CDS) of SlERF.D6 was obtained from the Sol Genomics Network database (see text footnote 1) and amplified from the cDNA of Micro-Tom leaves, then it was introduced into the entry vector pDONR207 and the destination vectors using the Gateway Cloning Technology (Curtis and Grossniklaus, 2003). Virus-induced gene silencing (VIGS) specific fragment (300 bp) was selected through the VIGS tool in the Sol Genomics Network database. After amplification and ligation to TRV2 vector (Orzaez et al., 2006), the plasmid TRV2-SlERF.D6 can be introduced into Agrobacterium tumefaciens strain GV3101 (Weidi, Shanghai, China).

Tomato fruits were harvested at Br5 and the fruit skin was ground into a fine powder in liquid nitrogen. Total RNAs were extracted by using a RNAprep Pure Plant Plus Kit (Polysaccharides & Polyphenolics-rich) (TIANGEN, Beijing, China). The total RNA (2 μg) was reverse-transcribed to cDNA using the TransScript^® One-Step RT-PCR SuperMix (TransGen, Beijing, China) according to the manufacturer’s instructions at 42°C for 30 min and 85°C for 10 s. A quantitative real-time PCR (qRT-PCR) was performed using the SYBR Premix ExTaq (Vazyme, Nanjing, China) by Quantstudio™ Real-Time PCR Software. Tomato ubiquitin gene (Solyc01g056940) was used for normalization, the relative expression level of each gene was calculated using the 2^–ΔΔCt method (Livak and Schmittgen, 2001). The primer pairs for qRT-PCR were designed using the Primer 5.0 and blasted at National Center for Biotechnology Information (NCBI) database to ensure primer specific. Three biological replicates were performed for each qRT-PCR analysis. Primers are listed in Supplementary Table 7.

For the measurement of SGAs, tomato fruits were harvested and half of the fruit pericarp was freeze-dried. Metabolite extraction and profiling were performed in the same way as described in the previous study (Zhu et al., 2018).

The plasmid TRV2-SlERF.D6, TRV2 and TRV1 were introduced into A. tumefaciens strain, then injected into MG fruit with the TRV1 vector combined with TRV2 or TRV2-SlERF.D6 plasmid, the TRV2 vector was used as control. The injection was performed as described previously (Orzaez et al., 2006; Fantini and Giuliano, 2016). The skin samples were collected 10 days after injection and immediately placed in liquid nitrogen. The vector pBin18-E8-SlERF.D6 used for transient overexpression was introduced into Agrobacterium and then injected into the tomato MG fruit in the same way, as the pBin18-E8 vector was used as control (Orzaez and Granell, 2009; Zhao and Zhao, 2009). Samples (skin) were taken into liquid nitrogen after 8 days.

The promoter fragments of GAME2, 5, 12, and 18 were amplified and introduced into pHis2 to construct bait vectors (BD) and the CDS of SlERF.D6 was introduced into pGADT7 to construct a prey vector (AD). All plasmids were transformed into Y187 Chemically Competent Cell (TOLOBIO, Shanghai, China) and cultured on SD/-Leu-Trp at 30°C for 2 days. The positive yeast strains were picked and diluted with sterile water to four concentrations. These suspensions were spotted on an SD/−Leu-Trp-His medium and incubated at 30°C for 3 days.

The CDS of SlERF.D6 was introduced into the pEAQ-HT-DEST2 vector acting as an effector vector and the promoter fragments of GAME12 were amplified and introduced into the modified pH2GW7 vector containing the firefly luciferase (LUC) gene and the Renilla luciferase (REN) gene acting as a reporter vector. All plasmids were transferred into A. tumefaciens strain GV3101 (Weidi, Shanghai, China) and then injected into the Nicotiana benthamiana leaves. The activities of LUC and REN luciferase were measured using the Dual-Luciferase Assay Kit (Promega) according to the instruction manual. The relative reporter gene expression levels were expressed as the ratio of firefly LUC to REN luciferase (LUC/REN). Six biological replicates for each sample were performed.

By conducting an integrative analysis on the tomato ITAG3.2 genome and multiple other databases, we identified a total of 176 AP2/ERF genes in tomato. Multiple protein sequence alignment and phylogenetic analysis revealed that all the AP2/ERF proteins in A. thaliana (127) and tomato (176) are mostly divided into three major clades (Supplementary Figure 1). Based on the above phylogenetic tree and the number of AP2 domains, 146, 26, and 3 of the genes that were identified as belonging to the ERF, AP2, and RAV subfamilies, respectively. In addition, a single gene, Solyc09g059510, that was homologous to At4g13040 (Giri et al., 2014) belonged to the Soloist subfamily. We explored the size of sequence, predicted subcellular localization, molecular weight, pI, and intron number of SlAP2/ERFs (Supplementary Tables 1, 2). In summary, new SlERFs were identified and their characteristics and structures herein were statistically determined.

Tomato breeding has undergone the two stages of domestication (PIM-CER, i.e., Solanum pimpinellifolium to S. lycopersicum var. cerasiforme) and improvement (CER-BIG, i.e., S. lycopersicum var. cerasiforme to S. lycopersicum). Fruit weight is artificially selected upon in modern tomato breeding programs, possibly leading to the selection of SlAP2/ERF genes (Lin et al., 2014). To characterize the role of these genes in domestication or improvement, we integrated the identification results of domestication/improvement sweeps in our previous study and the distribution of SlAP2/ERF genes (Zhu et al., 2018). If genes were located in the domestication or improvement sweeps, it means that these genes may regulate certain traits, such as fruit size that were selected by ancient farmers. We call these genes “domesticated” or “improved” genes and marked all domesticated or improved SlAP2/ERF genes on chromosomes (Figure 1). The results showed that one RAV gene and eight (more than 30%) AP2 genes were selected by humans. In total, 36 genes in the ERF subfamily located in the selection sweep, such as Solyc01g090340, Solyc10g009110, and Solyc03g119580 which have been reported to regulate the SGA synthesis, fruit growth, and development (Pirrello et al., 2012; Shi et al., 2012; Thagun et al., 2016). These data indicate that the SlAP2/ERF genes play a crucial role in tomato development and were introduced during the selection for fruit weight.

To further characterize the role of SlERFs in tomato, a phylogenetic tree was constructed with A. thaliana ERF proteins using the NJ method. These genes were divided into two major branches with nine groups. They were named DREB and ERF (Supplementary Figure 3), it is similar to other plants. Then, previously published tomato RNA-seq data collected from developing tomato fruit were used to explore SlERF-related expression patterns (Li et al., 2020). A total of 93 ERF genes with high transcription level (TPM > 1 at least one stage) were expressed from 10 days post-anthesis (DPA) to Br15 (days post-breaker stage) (Figure 2A). Hierarchical clustering showed that the ERF genes expression was divided into six clusters. During early fruit development, these ERF genes were highly expressed in cluster 4. During the late fruit development, they were more abundant in clusters 2, 3, and 6 (Figure 2B). These results imply that highly expressed ERFs play a stronger regulatory role during late fruit development.

In addition, we performed pairwise comparisons on RNA-seq data generated in this study during four key fruit development periods (IMG, MG, Br, and Ripe) to identify the DEGs. Among them, MG possessed the largest number of DEGs (4,803) compared with Ripe and 1,660 DEGs were identified between IMG and Br, which is the minimum number (Supplementary Table 4 and Figure 5). We found 14 genes that were not only significantly differentially expressed during tomato development, but also in an expression cluster we speculated was associated with fruit ripening (Figure 2B). Groups III and VII contained 28.6 and 21.5% of the 14 genes, respectively, which is more than any other group (Supplementary Figure 6). Several of these ERF genes, such as CBF1 and ERF1a, reportedly boosted the fruit development (Huang et al., 2004). These results suggest that our approach is appropriate and that some ERF genes with currently unknown functions may greatly influence the fruit development and ripening.

To understand whether other ERF genes can regulate SGA biosynthesis during tomato fruit development, the previously published tomato RNA-seq data were analyzed using WGCNA (Li et al., 2020; Supplementary Figures 7–9 and Supplementary Table 5). A total of 31 modules (each module is labeled with a unique color) were produced after filtering out low-transcription level (FPKM < 1) genes. Next, we associated each of the co-expression modules with metabolite concentrations in developing tomato fruits reported by previous studies (Li et al., 2020). Four of the modules were significantly positively correlated with SGA accumulation (PCC > = 0.60, p < 0.05), and three modules were significantly negatively correlated with SGA accumulation during tomato fruit development. In one of the modules, the gene expression was both positively and negatively correlated with SGA expression (Figure 2C). This suggests that the genes in these modules regulate SGA biosynthesis. Among them, it was found that the yellow and black modules were significantly positively correlated with Tomatidine-O-rhamnoside expression. Furthermore, we constructed the co-expression network between SlAP2/ERF and SGAs pathway-related genes using RNA-seq data by this study. When we looked for SlAP2/ERF and SGA pathway genes, we found 10 GAME genes, 5 sterol candidate genes, 28 upstream of genes involved in these pathways, and 29 SlAP2/ERF genes (PCC > 0.8 or PCC < −0.8 and p < 0.05). Furthermore, four SlAP2/ERF genes (SlAP2-004, SlRAV-002, SlERF.H21, and SlERF.D6) were highly co-expressed with multiple (gene number >10) SGA pathway genes (Figure 2D). This suggests that these four SlAP2/ERF genes regulate the expression of other SGA pathway genes which, in turn, regulates the SGAs biosynthesis. Since these four genes were also present in WGCNA modules that were highly associated with SGA pathways, we further speculated that they play important roles in the SGA biosynthesis.

In plant breeding, genetic gain depends on broad genetic variation (Shi and Lai, 2015). To find the genetic sources of natural variations in SGAs levels during tomato development, metabolite profiling data and high-quality SNPs were obtained from previous studies of 349 S. lycopersicum accessions (Supplementary Table 6). To identify the potential genetic basis of natural variations in SGA levels in tomato, mGWAS was performed on accessions displaying substantial differences in SGAs biosynthesis. One ERF gene locus, Solyc04g071770 (SlERF.D6), was significantly [p = 4.19 × 10^–7, Mixed Linear Model (MLM)] associated with the accumulation of esculeoside A (Figure 3A). Of note, this metabolite is produced downstream of the SGA biosynthesis pathway.

On the Manhattan plot, the most significant SNP was found at Chr.4:58746511 (Figure 3B). The LD of the SNPs in this region ranged from 57.75 to 59.875 Mb, as shown on the LD heatmap (Figure 3C). This result shows that there is an LD block near the most significant SNP, and SlERF.D6 (inferred to be a regulator of SGA biosynthesis) is in this block. These results indicate that SlERF genes may be responsible for the differences in SGA biosynthesis among the various tomato genotypes.

Further investigation revealed that the lead SNP, Chr.4:58746511, is a non-synonymous T to A mutation found within the coding region of Solyc04g071770. The SNP, which is located at the 593 bp position of the gene, converts an amino acid from Val to Glu. Haplotype analysis revealed that the SNP is significantly highly associated with the content of esculeoside A with metabolite ID of SlFM1998 (Figure 3D and Supplementary Table 6). Thus, it could be responsible for the variation observed in the developing tomato fruit. We found that the MAF of this SNP was not uniform among the tomato subpopulations. The frequency of the alternate allele was less than 3% in BIG, about 14.5% in CER, and as high as 54.5% in PIM (Figure 3E), suggesting that this gene was inadvertently selected by ancient farmers. To further understand their selection process, we measured the nucleotide diversity in this region of the genome (Figure 3F). Although it did not reside with a domestication or improvement sweep, a more pronounced difference between PIM, CER, and BIG genomes was found, hinting that it could be a target for future domestication sweeps. Based on these results, we found a natural variation of a gene regulating SGA biosynthesis, and the frequency of SNP was significantly different from wild tomato to cultivars.

To further understand the function of SlERF.D6, a qRT-PCR was performed to characterize the expression pattern of SlERF.D6 in different tomato tissues and different stages of fruit ripening. The result showed that the expression level gradually increased as fruit ripening progressed, indicating that SlERF.D6 may influence fruit ripen (Figure 4A). To investigate the role of SlERF.D6 in tomato fruit ripening, we performed transient expression in tomato fruit. In recent years, transient expression techniques have been applied to several species of plants (Hamilton et al., 1992; Liu et al., 2010). We first transiently overexpressed SlERF.D6 under the control of a fruit-specific E8 promoter (as shown in Experimental procedures). The qRT-PCR analysis showed that the expression level of SlERF.D6 was significantly higher in E8: SlERF.D6 tomato fruits than in the control sample (E8-EV). Thus, we successfully obtained three separate SlERF.D6 transient overexpression lines (Figure 4B and Supplementary Table 8). After injection, we observed that the fruit from the E8: SlERF.D6 lines split slightly more than the fruit from the control sample after 5 days. However, the color of fruit from the E8: SlERF.D6 lines were much redder relative to the control sample after 8 days (Figure 4C).

In addition, we performed VIGS to further explore the effect of SlERF.D6 on fruit development, which has been widely used in tomato fruit (Orzaez et al., 2006, 2009; Fantini and Giuliano, 2016). By TRV-mediated gene silencing, we obtained several lines with the reduced expression level of SlERF.D6 (TRV2: SlERF.D6) (Figure 4D and Supplementary Table 8). Besides, its phenotype was opposite to E8: SlERF.D6, when the fruit of EV (TRV2) had mostly turned red, TRV2: SlERF.D6 lines did not yet start to enter the stage of breaker (Figure 4E). These findings indicate that SlERF.D6 had a strong effect on the growth and development of tomato fruits, and that the overexpression of it promoted fruit ripening.

To better explore the relationship between SlERF.D6 and SGA metabolism, we measured the transcription level of several SGA biosynthesis-related genes using qRT-PCR analysis (Figures 5A–C). The results show that the expression levels of GAME2, GAME5, and GAME18 were higher in the E8: SlERF.D6 lines compared with the EV-E8 line, whereas the expression levels of GAME12 were lower (Figure 5B and Supplementary Table 8). In TRV2: SlERF.D6 lines, the expression levels of GAME2, GAME5, and GAME18 were lower and the expression levels of GAME12 were higher (Figure 5C). Additionally, we measured the concentration of SGA biosynthesis-related metabolites using liquid chromatography–mass spectrometry (LC-MS). The data showed that the levels of furostanol-26-aldehyde, β1-tomatine, α-tomatine, and esculeoside A were significantly higher in the E8: SlERF.D6 lines compared with the EV-E8 line, whereas the levels of tomatidine were lower (Figure 5D and Supplementary Table 9). Perhaps, the opposite was true in the TRV2: SlERF.D6 lines (Figure 5E). It is noteworthy that the level of esculeoside A was at least fivefold lower in the fruit from the TRV2: SlERF.D6 lines relative to the fruit from the EV-TRV2 line, which was significant than other metabolites (Figure 5E).

To investigate whether SlERF.D6 binds to the promoters of GAME genes to directly modulate SGA metabolism in tomato fruit ripening, a yeast one-hybrid (Y1H) assay was performed. The results showed that yeast containing the combination of SlERF.D6 and GAME12pro can grow normally on −Leu/Trp/His medium, and the control yeast contained the combination of pGADT7 empty vector and GAME12pro showed the opposite result. However, there was no similar result on the other three genes (Supplementary Figure 10A), which means that SlERF.D6 binds to the promoter of GAME12, but not to GAME2, GAME5, and GAME18. To further confirm the results, we performed a dual-luciferase experiment in tobacco leaves. The dual-luciferase analysis showed that co-expression with SlERF.D6 significantly decreased the luciferase activity of GAME12pro but had no effect on other GAME genes (Supplementary Figure 10B). These results suggest that SlERF.D6, acting as a transcriptional repressor, negatively regulates the synthesis of tomatidine via the direct transcriptional repression of GAME12 and it positively regulates the catabolism of tomatidine via the indirect transcriptional promotion of other three genes.