text:
Multi-Omics Analysis Reveals That SlERF.D6 Synergistically Regulates SGAs and Fruit Development

Steroidal glycoalkaloids (SGAs) are cholesterol-derived molecules that contribute to the pathogen defense in tomato but are toxic and considered to be antinutritional compounds to humans. APETALA2/Ethylene Responsive Factor (AP2/ERF) family transcription factors (TFs) play an indispensable role in various biological processes, such as plant growth and development, fruit ripening, biotic and abiotic stresses responses, and SGA biosynthesis. In this study, we identified 176 AP2/ERF genes that were domesticated or improved SlAP2/ERF in the tomato variome (Solanum lycopersicum) within either domestication or improvement sweeps, respectively. According to the RNA-sequencing data, 93 of the ERF genes with high transcriptional level (Transcripts Per Million, TPM > 1) belong to six clusters. Weighted gene co-expression network analysis (WGCNA) and metabolite-based genome-wide association study (mGWAS) analyses revealed that the expression level of the Solyc04g071770 (SlERF.D6) gene in the cluster six gradually increased as the fruit matured. Transient transformation verified that the overexpression of SlERF.D6 significantly promoted fruit ripening and regulated the expression of multiple genes in the SGA synthesis pathway, thereby affecting the SGA content of the fruit. Virus-induced gene silencing (VIGS) showed that the silencing of SlERF.D6 delayed fruit ripening and influenced the content of SGAs. Our data provide new insights into AP2/ERF TFs in tomato, offer a candidate TF for fruit development and steroidal glycoalkaloids, and provide new resources for tomato breeding and improvement.

INTRODUCTION
Tomato (Solanum lycopersicum) originated in the Andean region of South America. It is a model plant for genetic improvement and biological research system for Solanaceae crops, such as potato, pepper, and eggplant (Ranc et al., 2008). It is the most valuable fruit and vegetable crop in the world and contributes significantly to the amount of nutrients required to feed the human population (Klee and Tieman, 2018). Solanaceae plants produce a specific metabolite called steroidal glycoalkaloids (SGAs), the glycosylated forms of steroidal alkaloids (SAs), which help the plants protect themselves from bacteria, fungi, viruses, and certain insects during growth. For humans, SGAs are antinutritional compound that can either cause toxic reactions by inhibiting the cholinesterase activity in neurons or rupture lipid membranes through sterols (Aziz et al., 2012;Dzakovich et al., 2020;Hennessy et al., 2020).
In recent years, a series of research progresses have been made in the structural genes and transcription factors (TFs) of SGAs metabolic synthesis (Itkin et al., 2013). Glycoalkaloid metabolism 1, 17, 18, and 2 (GAME1, GAME17, GAME18, and GAME2) are glycosyltransferases, mainly performing multi-step glycosylation in the SGAs metabolic pathway. In addition to this, some P450s (GAME7, GAME8, GAME6, and GAME4), a dioxygenase (GAME11), and a transaminase protein (GAME12) are involved in the hydroxylation and oxidation of the cholesterol skeleton. Then, GAME5 is involved in the production of esculeoside A (Itkin et al., 2013). In addition, some TFs are involved in the regulation of SGAs. GAME9, an APETALA2/Ethylene Responsive Factor (AP2/ERF), plays a synergistic role in regulating the SGAs synthesis pathway. A 1 bp substitution in GAME9 was significantly associated with a high SGA content. It also affected the binding capacity of GAME17 promoter, a core SGA biosynthetic gene, with GAME9 (Cárdenas et al., 2016;Yu et al., 2020). This indicates that AP2/ERF plays a key role in the biosynthesis of SGAs.
In general, plant TFs play vital roles in plant growth, development, and environmental stress responses (Shu et al., 2016). AP2/ERF superfamily is one of the largest groups of TFs (Jofuku et al., 1994;Wessler, 2005), and the AP2/ERF proteins are involved in biological processes related to the plant growth and response to the biotic and abiotic stress in many species (Sharabi-Schwager et al., 2010;Qi et al., 2011). OsEATB, an AP2/ERF superfamily gene in Oryza sativa subsp., reduces rice plant height and panicle length at maturity by downregulating a gibberellin biosynthetic gene in the presence of ethylene (Qi et al., 2011). The DREB/CBF2 were shown to be involved in plant cold stress response, and the overexpression of CBF2 in Arabidopsis delayed plant growth, flowering, and senescence (Sharabi-Schwager et al., 2010). In tomato, SlAP2/ERFs have been shown to affect many biological functions, such as the regulation of growth of various phytohormones (i.e., cytokinin, ethylene, and methyl jasmonate) (Chung et al., 2010;Liu et al., 2013;Li et al., 2017;Sun et al., 2018;Lin et al., 2020;Zhang et al., 2021). Although Liu et al. (2016) updated the number of AP2/ERF genes to 146, there are still many potential AP2/ERF genes unknown. With the development of sequencing technology and multi-omics, it is crucial to reidentify AP2/ERF genes using new tomato genome data, and to analyze the function of AP2/ERF genes, such as the regulation of metabolism, in combination with multi-omics.
In this study, we re-identified the members of the SlAP2/ERF gene family using NGS tomato genome data, and the number was updated to 176. Additionally, we analyzed the characteristics and the evolutionary and phylogenetic relationships between these genes. By combining transcriptomic, metabolomic, and genomic, we had inferred and validated the functions of AP2/ERF genes, which revealed that some of them were selected during tomato breeding. Using both weighted gene co-expression network analysis (WGCNA) and metabolite-based genomewide association study (mGWAS), we found that SlERF.D6 coordinates the development of tomato fruit and the metabolic synthesis of SGAs. In this study, we identified the role of SlERF.D6 in fruit ripening and the synthesis of SGAs to provide a theoretical basis for cultivating better tomato varieties.

Plant Materials and Sampling
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 (RNAseq) 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.

Identification of Tomato AP2/ERF Gene Family
Tomato ITAG3.2 genome data were obtained from Sol Genomics Network. 1 The Hidden Markov model (HMM) data of AP2 domain were downloaded from Pfam database 2 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) 3 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. 4 The predicted subcellular localization data for all AP2/ERF genes were obtained via the GenScript website. 5

RNA Sequencing and Data Analysis
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 RNAseq 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 and Phylogenetic Tree Construction
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. 6 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.

Chromosomal Distribution, Gene Duplication, and Synteny Analysis
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. 7

Weighted Gene Co-expression Network Analysis and Co-expression Analysis
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.

Genome-Wide Association Analysis and Linkage Disequilibrium
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 genomewide 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 . 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 2 value by Plink (v1.9). 8

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

Expression Analysis by Quantitative Real-Time PCR
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 R 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 TM 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.

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

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

Yeast One-Hybrid Assay
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.

Dual-Luciferase Transient Expression Assay
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.

Identification of the SlAP2/ERF Genes
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.

SlERF Expression Regulates Fruit Ripening and Steroidal Glycoalkaloid Accumulation
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 . A total of 93 ERF genes with high transcription level (TPM > 1 at least one stage) were FIGURE 1 | Distribution of the identified 176 SlAP2/ERF genes across the tomato genome. The density of genes on 12 chromosomes are shown. All tomato chromosomes are drawn to scale based on their actual physical lengths. The genes with red color represent that they locate on domestication or improvement sweeps. The color of chromosomes represents the number of genes. expressed from 10 days post-anthesis (DPA) to Br15 (days postbreaker 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) FIGURE 2 | RNA-seq analysis of ERF genes during fruit development. (A,B) Expression pattern of the ERF genes during the fruit development of tomato. The gradient color (blue to red) in the heatmap represents the transcription level after normalization, and other colors represent different clusters. (C) Module-sample association. Each row corresponds to a module with different colors, labeled with a color as in Supplementary Figure 9. The color of each cell at the row-column intersection indicates the correlation coefficient between the module and the steroidal glycoalkaloids (SGAs) content, and the modules with black border represented that they are highly correlated with SGAs. (D) Analysis of co-expression between SlAP2/ERF genes and other SGA related genes. Different colored circles represent different genes, red is the AP2/ERF identified in this study, and different circle sizes are the number of their high correlation with other genes.
Frontiers in Plant Science | www.frontiersin.org 6 April 2022 | Volume 13 | Article 860577 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 ; 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 lowtranscription level (FPKM < 1) genes. Next, we associated each of the co-expression modules with metabolite concentrations in developing tomato fruits reported by previous studies . 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-Orhamnoside expression. Furthermore, we constructed the coexpression network between SlAP2/ERF and SGAs pathwayrelated 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.

Use of Multi-Omics to Identify Steroidal Glycoalkaloid Biosynthesis Candidate Genes
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.

SlERF.D6 Promotes the Ripening of Tomato Fruit
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., 2006Fantini 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.

SlERF.D6 Regulates Steroidal Glycoalkaloids Metabolism During Tomato Fruit Development
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 contents of some related metabolites in TRV2: SlERF.D6 lines. The Y -axis represents relative content. All the above error bars represent the SD (n = 3) (*p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001; Student's t-test). 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 coexpression 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.

DISCUSSION
Transcription factors, such as the NAC, MYB, bHLH, WRKY, and AP2/ERF families, are widely involved in the fruit development and metabolic regulation in plants (Groszmann et al., 2008;Machemer et al., 2011;Liu et al., 2018). In recent years, there are numerous reports on the identification and analysis of AP2/ERF gene families in plants, such as tomato, rice, rape, and other plants. They can regulate fruit ripening and respond to a variety of stressful environments (Du et al., 2016;Klay et al., 2018;Donde et al., 2019;Song et al., 2020;Shoji and Yuan, 2021). With the development of sequencing technology, more and more plant genome data are reassembled to obtain higher quality genome (Jung et al., 2019;Bohra et al., 2020), and the research of AP2/ERF family is more and more in-depth and comprehensive. Thus, identifying AP2/ERF TFs can considerably improve our understanding about their evolution and function in various plant species. A total of AP2/ERF genes in the tomato genome resulted in the identification of 176 members, and previous studies have revealed that with the AP2/ERF TFs, the number was 146 and we expanded the number of this family and given the basic characterization (Pirrello et al., 2012). In addition, we found 45 of them under domestication or improvement stage, it suggested that these genes may regulate some traits that were relevant to humans.
Tomatidine that determine the bitterness are gradually synthesized into esculeoside A accumulated in red mature fruits (Nakayasu et al., 2018). We performed mGWAS to scan the gene SlERF.D6 associated with esculeoside A (Figure 3A). Next, we found during fruit development, the content of tomatidine decreased rapidly, while esculeoside A content increased significantly (Figures 5D,E), which was consistent with previous studies (Cárdenas et al., 2016). In addition, the ERF.D6 presented in the group which might be the main factor for the change of fruit color and firmness, and genetically modified breeding will be generated by genetic engineering technology to improve the ripening time and nutritional quality of tomato in future. Therefore, this study conjectured that SlERF.D6 could regulate SGA metabolism and fruit development in addition to the functions of resisting abiotic stress .
Natural variation in ERF.D6 can disrupt the normal metabolic levels of SGAs, and we also found it under the domestication or improvement selection (Figure 3E), it showed that ERF.D6 underwent natural variation in the population and the variation in this gene was associated with esculeoside A, and this different trait was selected by ancient farmers. In addition, the natural variation in ERF.D6 can promote fruit ripening, allowing farmers to obtain ripe tomato fruit faster and shorten the growing cycle. Thus, the candidate gene ERF.D6 both regulates fruit development and at the same time is an important gene for altering the nutritional varieties of tomato (reducing toxicity). Selection for haplotype of ERF.D6 and rational crosses could lead to new tomato varieties with faster fruit ripening and less antinutritional metabolites, which could be of great use for tomato breeding applications and crop genetic improvement.
Plant growth and development are coordinated by a complex network of interacting phytohormones, and the interplay between ethylene and auxin signaling is critical for a wide range of plant developmental processes. In a previous study, Liu et al. (2018) demonstrated the role of tomato ERF.B3 in integrating ethylene and auxin signaling by directly regulating Sl-Aus/IAA27, while Wang et al. (2009) commented on the role of Aux/IAA in plant specialized metabolism. In another work, it showed that the reciprocal AP2/ERF gene clusters regulate the biosynthesis of specific metabolites in a variety of plants (Paul et al., 2020), such as ORCA5 overexpression in tobacco and NICOTINE2 ERF189 overexpression in Catharanthus roseus activate nicotine and terpenoid indole alkaloid pathway genes, respectively (Kato et al., 2014). In SGAs biosynthesis, studies shown that phytohormones, such as JA and GA, can regulate SGAs to resolve chemical defense with growth in tomato (Nakayasu et al., 2018;Panda et al., 2022). In view of the above research content, we speculate that SlERF.D6 may be induced by phytohormones to regulate SGAs.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: National Center for Biotechnology Information (NCBI) BioProject database under accession number PRJNA783378.

AUTHOR CONTRIBUTIONS
SW conceived the project and supervised the study. HG and YD conducted the data analysis and curation. MM and PC contributed to experimentation and methodology. RC and JL performed the metabolite profiling. YZ, LS, CW, CL, YL, QB, and TT participated in the material preparation. HG, MM, and PC wrote the first draft. SW and JY reviewed, edited, and finalized the manuscript. All authors discussed the results and commented on the manuscript.


provenance: AgriculturalAndFoodSciences,2022-2022-Agricultural And Food Sciences-2022.jsonl.gz.gz

text:
Classification of Sidr honey and detection of sugar adulteration using right angle fluorescence spectroscopy and chemometrics

Sidr honey is vulnerable to adulteration with low-grade honey and sugar syrups, which compromises its nutritional and medicinal value, demanding fast and reliable analytical tools for quality assessment. In this study, fluorescence spectroscopy was employed to assess the quality of a honey samples, specifically, Sidr, unifloral (Acacia) and multifloral (Acacia, Carisa and Justicia) honey. Fluorescence spectroscopy revealed characteristic spectral signatures of Sidr honey, compared to Acacia and multifloral honey. In addition, cane sugar syrup was artificially added to Sidr honey at different concentrations. These spectral signatures were exploited for the machine-assisted classification of Sidr, sugar syrup and different concentrations of Sidr–sugar mixture. The bagging classification algorithm generated values of sensitivity and specificity close to unity, indicating its ability for efficient discrimination of the samples. Fluorescence spectroscopy in tandem with chemometrics could potentially be used as a rapid analytical tool to identify Sidr honey and its sugar adulteration. Supplementary Information The online version contains supplementary material available at 10.1007/s00217-022-04008-9.

Introduction
Honey has been a part of human food since centuries as a nutrient as well as therapeutic agent. Nowadays, in the period of COVID-19, many people are consuming it in combination with lime or green tea due to its known antiinflammatory and antimicrobial effects. The honey consists of carbohydrates, mainly monosaccharides (fructose and glucose) and minute amount of disaccharide (sucrose) and other types of sugars (oligo and tetrasaccharides). Besides, it contains a small amount of protein, phenolic substances, HMF (hydroxyl methyl furfural), vitamins, minerals and water. The composition of honey varies depending upon floral origin (plant type), climatic variation, geographical region and contribution by beekeepers (Ankalam 1998). Honey from a defined botanical origin is unifloral honey and is distinctive from multifloral honey due to specific organoleptic properties [26]. The anti-inflammatory and antioxidant properties of honey are due to its phenolic compounds, flavonoids, carotenoid derivatives and organic acids [1].
Among all honey, Sidr honey is expensive and Yamen is famous for its Sidr honey since centuries due to its nutritional and therapeutic qualities. Sidr is also known as Jujuba, these tress grow in the desert areas of Yamen, belonging to the genus Ziziphus, are the major source of collecting nectar for pollinating bees. In Pakistan, seven species of Ziziphus have been reported and jujuba trees are grown mainly in Karak region. Sidr honey is collected from this region, from where it is exported to Arab countries. Pakistani flora is represented by Ziziphus species from Sindh, Pujab, Khyber Pakhtunkhwa (KPK) and Balochistan provinces [9].
According to the European Council directive (Codex 2001; Council Directive 2001/110/EC), honey is a natural sweet substance produced by honeybees from the nectar of plants and transformed by combining with specific substances of their own. Honey is adulterated directly or indirectly with inexpensive sugar syrups like high fructose corn syrup, corn syrup, inverted sugar syrup, sugarcane syrup and jaggery syrup. Direct adulteration is done by adding sugar syrup to honey after production and indirect adulteration is done by feeding sugar syrup to honeybees during nectar flow period to get more honey [11]. Commercial sugar is extracted mainly from sugar cane or sugar beet and is mainly composed of sucrose (99.8%) and 0.2% of polysaccharides, organic and inorganic non-sugar compounds and colorant polymers [5].
Different methods have been reported to detect sugar adulteration in honey [25,28] like chromatography (HPAEC-PAD) [8], thin-layer gas chromatography (GC-MS) [20], carbon isotope ratio (13C/12C) technique [11] and nuclear magnetic resonance (NMR) [17]. These classical methods are tedious, costly and destructive with several limitations [4]. Fluorescence spectroscopy being non-destructive and rapid food analytical technique is known for its wide use in food characterization and determining molecular components, additives and adulterants [21]. Intrinsic fluorophores in food absorbs ultraviolet or visible light at a specific wavelength and emit fluorescence at higher wavelength [14].
Chemical profiling, quality assessment and botanical origin of honey have been reported using fluorescence spectroscopy in combination with multivariate analysis. Front face or reflectance fluorescence spectroscopy (FFFS) has been used by many researchers due to the viscous and opaque nature of honey [13,18,19]. Classical right angle fluorescence methodology was also evaluated to discriminate honey samples solubilized in buffer and methanolic solution from the region of Croatia [24]. In these studies, fluorescence spectroscopy coupled with PCA and different regression models were applied. Owing to the complex nature of honey, different statistical tools for the characterization of honey and adulterant in it are required. Normally, unsupervised multivariate data analysis tools such as PCA or hierarchical cluster analysis (HCA) are used. In this study, in addition to PCA, we performed ensemble learning to achieve high accuracy for classification of dataset. An ensemble method combines the predictions from multiple machine learning models for more accurate predictions as compared to any individual model. Bootstrap aggregation (or Bagging for short) is a very powerful ensemble method that independently trains several machine learning algorithms using randomly chosen subsets of data. It reduces the variance for algorithms with high variance [16].
The research in our group was started with the aim to establish quick and non-destructive analytical test for characterization of Sidr honey and detection of adulteration in it. Sidr honey collected from various beekeepers and villagers was fingerprinted and distinctive fluorescence spectral signatures were identified in comparison with other monofloral honey samples [2] using synchronous fluorescence spectroscopy (SFS) and PCA. SFS is the multicomponent fluorophore analysis as compared to right angle fluorescence spectroscopy [10].
In the present work, application of right angle fluorescence spectroscopy along with machine learning method has been applied to identify a single wavelength that can be used quickly to fingerprint Sidr honey and sugar adulteration in it. Development of quick, eco-friendly and non-destructive analytical tool may also be useful for the development of handheld device for honey researchers, producers, exporters, consumers and retailers.

Sample preparation
The study was conducted on 25 samples of Sidr honey (S1-S25) along with 25 samples of Acacia (A1-A25) and 25 samples of multifloral (Acacia, Carisa and Justicia) honey (MF1-MF25), provided by Honey Bee Research Institute, Pakistan Agricultural Research Council, (PARC), Islamabad, Pakistan. Physiochemical, microscopic and sensory analysis of these honey samples was performed at PARC and listed in Supplementary Tables 1, 2, 3. The adulterated set included ten solution mixtures. Sugar syrup was prepared by adding 50 g of cane sugar in 50 mL of distilled water and heated in a water bath at 50 °C to dissolve sugar crystals. Different concentrations (5,10,15,20,25,30,35,40, 45 and 50%) in weight were prepared by mixing sugar syrup in pure Sidr honey. The Sidr-sugar mixtures were then heated in a water bath and vortexed for homogenous mixing.

Fluorescence spectroscopy
Spectra of all different types of honey including Sidr, Acacia, multifloral and Sidr-sugar mixtures were analyzed by fluorescence spectroscopy using FluoroMax®-4 (HORIBA Scientific, Jobin Yvon, USA) spectrophotometer fitted with a xenon lamp and a Fluorescence TM software. Right angle geometry was used for samples in a 1 cm quartz cuvette for recording emission spectra ranging from 365 to 690 nm. Two excitation wavelengths, i.e., 350 and 405 nm, were used with excitation and emission slits set at 2 nm.

Data processing and chemometrics
The spectra were processed and vector normalized using Matlab software (2014a). The fluorescence data of honey samples was subjected to principal component analysis (PCA) to differentiate groups using Origin software (Orig-inPro 2017, OriginLab.Corp., Northampton, MA).
In this study, the bootstrap aggregating-bagging for short-algorithm, was used for machine-assisted classification of Sidr, sugar and mixed sidr-sugar samples. A total of seven different classes were considered-one each for Sidr and sugar and five mixtures of Sidr-sugar. The complete cohort of data comprised 140 samples (i.e., 20 for each of the seven classes). The predictive capabilities of the bagging technique are leveraged on the famous classification and regression tree (CART) algorithm, where the input data space is transformed into a classic binary decision tree, as detailed below: • The position of maximum intensity signal in the fluorescence spectra (i.e., max ) collected for each sample was used as input data for the automated multiclass classification. • The data elements from the original input data space D were transferred to a new, equal dimensional data space D n ; this transformation was a random process so that some data points were not transformed at all while others were transformed repeatedly. • Based on the characteristics of the dataset D n , the algorithm generated multiple classifiers H n , where n = 1, 2,…, N (= number of characteristics of D n ). • Training followed by testing of each H n was executed on the dataset D n so that H n ∶ D n → R.

Results and discussion
The normalized fluorescence emission spectra of all investigated honey samples (Sidr, Acacia and mutifloral) with both excitation wavelengths (350 and 405 nm) are presented in Figs. 1, 2. The emission from pure Sidr honey can be seen in the range of 494-508 nm with λ exc 350 nm, whereas similar spectra was observed for Sidr with λ exc 405 nm (Fig. 2). On the other hand, MF and Acacia honey exhibited overlapping Emission maxima shifted to longer wavelengths from 474 to 490 nm in both Acacia and MF honey when excitation wavelength was changed from 350 to 405 nm. Such shifts to longer wavelengths point to emission from increased molecular size of fluorophore [12]. However the fluorescent cluster in pure Sidr honey remained unchanged at both excitation wavelengths.
The fluorescence spectral profile in honey is due to the contributions of main fluorophores, such as phenolic acids, flavonoids, amino acids, and vitamins. Thus, the presence of slightly broader hump in MF and Acacia as compared to Sidr honey represents overlapping fluorescent compounds owing to multifloral nature of sample predominantly from Acacia. Acacia plant was one of the honey sources in MF samples and seems like all the three sources (Acacia, Carisa and Justicia) shared the same fluorophores and thus showed overlapping emission spectra. Research studies have confirmed the presence of polyphenols in honey that are responsible for sensory properties (color, taste and odor) as well as antioxidative, antimicrobial, antitumoral and anti-inflammatory properties [15] depending upon the floral origin from which they are extracted.
Previously, we found a blue green fluorescence emission represented by broad hump in the range of 395-459 nm in commercial multifloral samples and 420-445 nm from samples provided by National Agriculture Research Council farms linked with HBRI using synchronous fluorescence spectroscopy (λ exc range 250-450 nm, Δλ = 60 nm). In Sidr honey, the main difference from the multifloral honey can be seen in the spectral region at 397 and 451 nm, represented by phenolic acids [2]. Synchronous fluorescence spectra reduce spectral overlaps and thus display broad fluorescence bands that contain chemical, physical and structural information of components. Alzubier & Okechukwu [3] reported that Sidr honey from Yemen has different phytochemicals such as flavonoids, steroid, alkaloids, saponins and tannins and possess strong anti-inflammatory, analgesic and antipyretic effects.
Though it was not difficult to compare the fluorescence emission spectra by visual inspection, still data were subjected to chemometrics for authentication. Principal component analysis is applied to reduce the dimensionality of data and analyze statistical differences between spectral data of different samples [23]. The resulting score plot for PCA applied on the normalized data of honey samples is presented in Fig. 3. The total variance for the two components, i.e., PC1 and PC2, corresponds to 93.4% and 5%, respectively. Sidr honey formed a cluster that is distant from Acacia and MF samples, suggested due to the presence of different phenolic compounds represented by longer emission wavelength. Most of the samples from Acacia and MF showed overlapping that can be explained by the presence of similar fluorophores in these honey types.
Sidr honey mixed with different concentrations (5,10,15,20,25,30,35,40, 45 & 50%) of sugar syrup was analyzed using excitation wavelength of 350 nm. We used 350 nm instead of 405 nm, as more distinct emission spectra were observed relating to Sidr and non-Sidr (MF and Acacia) honey types. Figure 4 shows the specific qualitative display of Sidr honey fluorescence emission maxima at 507 nm with gradual shift toward shorter wavelength subsequent to the addition of sugar syrup in honey. In this study, sugar syrup was expressed with a fluorescence emission band with maxima at 435 nm, also reported by Ghosh et al. [10] with excitation at 340 nm and emission at 430 nm that was suggested to be generated from nicotinamide adenine dinucleotide (NADH). Baunsgaard et al. (2006) analyzed beet sugar and cane syrup just before sugar crystallization and found emission from amino acid and colored precursors with higher wavelength in cane sugar as compared to beet sugar. The fluorescence in sugar syrup in visible wavelength might be generated from colored compounds (Maillard reaction products) during thermal processing of cane sugar [6] (Fig. 4).
PCA was applied to the five selected sugar adulterated concentrations of Sidr honey that showed discrimination of pure honey subjected to adulteration with cane sugar syrup. The score plot (Fig. 5) defined by two components is represented by a separate cluster consisting of Sidr and adulterated samples (10-30% sugar syrup adulteration) with dominance of PC1 (89.7%). The pure honey samples were well separated from adulterated samples (40% and 50% sugar syrup adulteration) and pure sugar syrup at PC1, but there was small variation for separation of adulterated samples up to 30%. The remaining component of the twodimensional score plot, i.e., PC2, showed small contribution A discriminant analysis based on bagging algorithm was also developed to differentiate pure and adulterated Sidr honey. The ability of the bagging classifier to handle multifaceted data, while concurrently diminishing the noise, variance and bias in the classification, has motivated its widespread use in applications pertaining to biomedical optics [22,27]. Moreover, the implementation of the bagging algorithm for data containing gaps, i.e., missing data points, is equally efficient in terms of accuracy. In perspective of these strengths, a custom version of the bagging algorithm was utilized in this study. Table 1 presents a summary of the classification of the samples and the quantitative metrics used for performance evaluation of the bagging classifier. Specifically, the algorithm was modeled and executed for classification of sevenclass problem-one class each for Sidr and sugar and five classes for different Sidr-sugar mixtures. which is subsequently exploited to compute different metrics that identify and quantify the classifier's efficiency. In particular, the range of sensitivity and specificity values of the bagging classifier for the seven classes was 0.55-1.0 and 0.94-1.0, respectively. For pure Sidr and pure sugar, the bagging classifier generated the highest values for these two metrics (i.e., sensitivity and specificity), indicating that the machine-based differentiation of samples in these two classes is substantially accurate. A relatively low values of sensitivity (= 0.55) and specificity (= 0.95) was observed for sugar 40. Other metrics computed for the appraisal of the classier included PPV, NPV, F-measure and MCC; the values of these metrics for all seven classes have been given in Table 1. Albeit all metrics demonstrated efficient performance of the bagging classifier (as the values are close to unity) for the seven-class problem, the details of specific trends related to each individual metric can be analyzed from Table 1. In future, such model will also be developed for other adulterants in Sidr honey such as fructose corn syrup, beet sugar syrup and glucose corn syrup.

Conclusion
The study revealed that right angle fluorescence spectroscopy, being a non-destructive and rapid technique, is a promising tool to fingerprint Sidr honey with no sample processing. This fast and simple technique along with chemometrics can be practically applied for the identification of Sidr honey at massive scale and detection of direct sugar adulteration in it. Moreover, this technique can be utilized in the development of handheld device that will be very beneficial not only for honey researchers, but also for producers, exporters, consumers and retailers for the authentication of Sidr honey.