The nucleotide and amino acid sequences of the AtWRKY and AtbHLH proteins of Arabidopsis were downloaded from The Arabidopsis Information Resource (TAIR, https://www.arabidopsis.org/). Nucleotide sequences of the SlWRKY and SlbHLH proteins of S. lycopersicum were retrieved from Phytozome (https://phytozome-next.jgi.doe.gov/) and Sol Genomics Network (SGN, https://solgenomics.net/), while their corresponding protein sequences were downloaded from the Plant Transcription Factors Database (PlantTFDB, http://planttfdb.gao-lab.org/). The subcellular localizations of SlWRKY and SlbHLH proteins were predicted using CELLO (http://cello.life.nctu.edu.tw/) ( Table S1 ). The protein sequences of putative SlWRKY and SlbHLH members were subjected and analyzed with Expasy ProtParam (https://web.expasy.org/protparam) and SGN (https://solgenomics.net/) to determine their theoretical isoelectric point (pI), molecular weight (Mw), and number of amino acids (AAs).

Neighbor-joining (NJ) phylogenetic trees of the 72 Arabidopsis with 81 S. lycopersicum WRKY and 153 Arabidopsis with 161 S. lycopersicum bHLH proteins were constructed using the software MEGA 11 with 1,000 replicates for bootstrap (BS) analysis for statistical reliability to compare the evolutionary relationships of WRKY and bHLH TFs across the species. We further performed NJ phylogenic tree analysis among WRKY and bHLH protein sequences in S. lycopersicum to test the reliability of the results.

The CDS and genomic sequences of the WRKY and bHLH of Arabidopsis and S. lycopersicum were retrieved from TAIR (https://www.arabidopsis.org) and Phytozome (https://phytozome-next.jgi.doe.gov/), respectively. The Gene Structure Display Server (GSDS, http://gsds.gao-lab.org/) web tool was used to analyze the structure of the members of the WRKY and bHLH gene families and to detect their exon–intron organization, by aligning the CDS sequences with genomic sequences. To identify the conserved motifs in these protein sequences, the Multiple EM for Motif Elicitation (MEME, https://meme-suite.org/meme/tools/meme) online server was used with the following parameters: number of repetitions, any; maximum number of motifs, 15; and optimum motif width set to ≥6 and ≤200 amino acid residues.

In order to determine the physical locations of the SlWRKY and SlbHLH genes in the S. lycopersicum genome, the starting and ending positions of all the identified genes on each chromosome were obtained from the Solanaceae family database (Sol Genomics Network). The MapInspect software (http://mapinspect.software.informer.com/) was used to map the physical locations of the SlWRKY and SlbHLH genes on their respective chromosomes. The duplicate chromosomal blocks were retrieved using the Plant Genome Duplication Database (PGDD, accessible at http://chibba.agtec.uga.edu/duplication/), and the WRKY and bHLH genes within the duplication block were identified. This allowed us to find duplicate S. lycopersicum WRKY and bHLH genes (Lee et al., 2012). Genes separated by five or fewer gene loci in a range of 100 kb distance were considered to be tandem duplicates, and those which were co-paralogs and located within duplicated chromosomal blocks in multiple locations as a result of duplication and chromosome rearrangement and shared >90% sequence identity were considered as segmental duplicates (Wei et al., 2007). The PGDD, a public web service database, was used to identify and characterize the genes in terms of intra- or interplant genomic syntenic relationships.

In the current study, tomato seeds (S. lycopersicum cv. Yegwang) were used. First, the seeds were surface-sterilized with 10% hypochlorous acid and 70% ethanol and washed with autoclaved distilled water to remove the impurities. The soaked seeds were germinated on hygiene filter paper in an incubator at 30°C in dark conditions. After successful sprouting, the seeds were planted in plastic pots in a greenhouse. The greenhouse temperature was kept constant at 28°C ± 2°C, 55% ± 5% relative humidity, and 16-h light/8-h dark photoperiod. The experiment was conducted utilizing three groups of plants: a) distilled water-treated plants (control plants), b) 1-mM Cd-treated plants, and c) 2-mM Cd-treated plants. After every 3 days, the plants were treated with their respective treatments under the same growth conditions. Four replicates were prepared per treatment. Leaf samples were collected randomly and immediately placed in liquid nitrogen and stored at –80°C in a fridge until analysis.

RNA was isolated using the RNeasy Plant Mini Kit (QIAGEN, Hilden, Germany). The quantity of total RNA was adjusted to 10 μg, and the NanoDrop ND-1000 (Thermo Scientific, Waltham, MA, USA) was used to evaluate the quality. Sequencing and analysis of libraries from three different biological replicates of each treatment were performed. The Illumina HiSeq 2000 platform was utilized to perform library construction according to the described approach (Liao et al., 2015), resulting in single-end reads with 51 bp. To identify variations in gene expression regulation between the cadmium-treated and non-treated plants, a computational pipeline of optimized tools was employed. For trimming and quality check, Trim Galore (Krueger, 2012) and FastQC (Andrews, 2010) were used, respectively. To align the reads to the reference genome, HISAT2 (Sirén et al., 2014) was used before the read count quantification by using Feature Count (subread_v2.0.2). The R software using the DESeq2 (Love et al., 2014) package was used for differential gene expression analysis.

The STRING protein interaction database version 11.5 (https://string-db.org) was used to identify the interacting protein networks and functional annotations.

The online software agriGO analysis toolbox (http://bioinfo.cau.edu.cn/agriGO/) was used to enrich gene ontology (GO) categories (Du et al., 2010) using the TopGO “elim” algorithm (Alexa et al., 2006) covering the following characteristics: biological processes, molecular processes, and cellular processes.

For qRT-PCR analysis, approximately 10 highly expressed genes (5 WRKY and 5 bHLH) from RNA-seq data were selected to authenticate the RNA-seq results. The Primer3 (https://bioinfo.ut.ee/primer3-0.4.0/) program was used to design primers for each selected gene as listed in Table S3 . The standard cDNA was produced using PCR Biosystems’ qPCRBIO cDNA Synthesis Kits after total RNA was diluted to a final concentration of 100 ng/μl, and transcript quantification was done as reported previously (Asaf et al., 2022). Actin (a housekeeping gene) was used as an internal control to normalize gene expression, and the comparative ΔΔCt method of qRT-PCR was utilized to calculate the expression level of the genes in control plants in comparison with Cd-treated ones.

In order to systematically identify and analyze WRKY and bHLH genes in the tomato genome, their sequences were downloaded as described above. Finally, a total of 81 and 161 non-redundant putative WRKY and bHLH genes were confirmed by their specific domain using SMART (http://smart.embl-heidelberg.de) ( Table S1 ). Subsequently, the downloaded gene sequences were compared with their expressed sequence tag (EST) sequences in PlantTFDB (http://planttfdb.gao-lab.org). Arabidopsis WRKY and bHLH genes were obtained from the TAIR (https://www.arabidopsis.org) and PlantTFDB (http://planttfdb.gao-lab.org). A total of 72 WRKY and 161 bHLH annotated TFs were extracted from these sources ( Table S1 ).

In order to determine the evolutionary and phylogenetic relationship among the S. lycopersicum WRKY as well as among bHLH proteins, unrooted neighbor-joining phylogenetic trees were built by a multiple sequence alignment of their amino acid sequences using the ClustalW program. As shown in Figure 1A , WRKY TFs of S. lycopersicum were divided into six main groups (group I to group VI). Notably, the maximum number, i.e., 19 (23.4%) of the WRKY members, was found in group I, followed by groups III and IV with 18 (22.2%) members. The least number, i.e., 6 (7.4%) of the WRKY members, was found in group V. Similarly, bHLH TFs were divided into eight main groups (group I to group VIII) based on the phylogenetic tree. The maximum number, i.e., 32 (19.9%) of the bHLH members, was found in group I, followed by group VIII with 30 (18.6%) and group V with 26 (16%) members. The least number, i.e., 7 (4%) of the bHLH members, was found in group II ( Figure 1B ). Furthermore, most of the members from the same phylogenetic group have the same number of exons and conserved exon–intron structure ( Figures 2 , 3 ). For instance, in groups III and IV, 18 (75%) and 14 (78%) of the SlWRKY members have three exons, respectively. On the other hand, in the phylogenetic tree of SlbHLH TFs, group VIII contained the majority of the genes with 6–13 exons. Our results also showed that most of the evolutionarily related members have similar motif compositions and the same subcellular locations ( Figures S1 , S2 and Table S1 ). The similarity in these features may be related to their specific physiological functions in tomato cells. The high bootstrap values (95%–100%) indicate a strong phylogenetic relationship between some pairs of the SlWRKY TFs. For instance, Solyc05g050050–Solyc05g050060, Solyc01g058540–Solyc04g051540, and Solyc05g050330–Solyc05g050340 pairs of SlWRKY genes showed a paralogous relationship, supported by 99% of bootstrap values ( Figure S3 ), while some pairs of SlbHLH TFs, i.e., Solyc01g014910–Solyc09g018130, Solyc01g020170–Solyc01g107970, Solyc03g119390–Solyc12g036470, Solyc09g089870–Solyc10g006510, and Solyc10g006510–Solyc08g075090, had 100% bootstrap values as shown in Figure S4 .

To identify the evolutionary history and phylogenetic relationship between the SlWRKY (81 members) and AtWRKY (72 members) and between the SlbHLH (161 members) and AtbHLH (153 members), unrooted NJ phylogenetic trees were derived from the alignment of the amino acid sequences of these TFs. The tree inferred from the WRKY TFs of Arabidopsis and S. lycopersicum organized WRKY into nine main groups (group I to group IX) ( Figure 3A ). Among the nine main groups, group IX is the largest one containing 29 members, followed by group I with 24 members and group IV and group VIII with 20 members, while group VI was the smallest which contained only nine members. For bHLH TFs, 11 groups (group I to group XI) were defined based on their constructed phylogenetic tree. Group V is the largest one with 50 bHLH members, followed by group I with 45 members and group III with 41 members, while group II has the least number (only 10) of members ( Figure 3B ).

Genes that descended from a common ancestral gene are likely to have the same functions. Therefore, homological analysis is widely used to predict the function of an uncharacterized gene (Dai et al., 2020). Paralogs are homologous genes evolved by duplication of an ancestral gene which can keep the same or diverge to a different functional role (Gabaldón and Koonin, 2013), and orthologous relationships between genes from different species can reflect the same evolutionary and functional niche in different species (Fang et al., 2010). The tree topology and the organization of most of the groups and subgroups resembled those from the S. lycopersicum and Arabidopsis individual trees. The trees presented in Figures 1 , 3 identified putative paralogs and orthologs. A total of 27 and 50 orthologous pairs of WRKY and bHLH genes, respectively, were identified across the species. Most S. lycopersicum and Arabidopsis paralogous genes observed in Figures 3A, B were already displayed as paralogs in the respective trees ( Figures 1A, B ), for instance, Solyc01g058540–Solyc04g051540, Solyc05g050330–Solyc05g050340, and Solyc05g050050–Solyc05g050060 pairs of the SlWRKY TFs and Solyc09g089870–Solyc10g006510, Solyc08g075090–Solyc08g075110, and Solyc01g020170–Solyc01g107970 pairs of the SlbHLH TFs. The exception of a few TFs for not displaying in the same paralogous pairing in both individual and combined trees may be due to the existence of an apparent ortholog from a different species, providing better support for the new location, e.g., Solyc12g014610–AT1G13960 and Solyc02g021680–AT2G40740 pairs of group VI of the combined phylogenetic tree of SlWRKY and AtWRKY TFs. It is worth noting that most of the described phylogenetic groups and subgroups were supported by additional evidence, such as gene structure, conserved motifs, gene ontology, protein–protein interaction, and similarity in physiological functions of most of the characterized genes.

Analysis of the orthologous and paralogous relationships is widely used to predict the function of an uncharacterized gene because the genes that descended from a common ancestral gene are likely to have the same functions (Dai et al., 2020). When comparing multigene families of different species, it is frequent and interesting to find multiple genes in one species that are orthologous to a single gene in the other, indicating recent duplications specific to the former. The orthologous genes with similar sequences and the same expression patterns in various organisms suggest a possible case of functional redundancy. There are no detailed data on the specific physiological roles of most S. lycopersicum genes. Since both WRKY and bHLH TFs have important roles in several important major developmental and physiological processes and responses to different environmental stimuli, it is critical to conduct extensive research on the WRKY and bHLH gene families in S. lycopersicum. The comparative genome-wide study of S. lycopersicum WRKY and bHLH with WRKY and bHLH TFs of Arabidopsis, which is the most extensively studied plant, will allow us to predict their physiological functions. For instance, AT1G01260 is involved in negatively regulating jasmonic acid during leaf senescence and making an orthologous pair by a strong BS score with the Solyc01g096050 gene of the bHLH family. Similarly, AT1G26945 is involved in ABA and salt stress responses resulting from an orthologous pair by a strong BS score with Solyc06g070910 of the bHLH family (Hao et al., 2021). So, based on these orthologous relationships, we can predict the function of Solyc01g096050 and Solyc06g070910 genes in S. lycopersicum.

To further understand the evolutionary patterns and gene duplication events, structural and compositional analyses of genes can be used as supporting evidence (Khan et al., 2021a). The WRKY and bHLH genes’ exon–intron patterns and conserved motifs were studied to provide insight into the evolution of these two gene families in the S. lycopersicum genome. The results showed that the number of exons and introns in WRKY genes ranged from 6 to 1 and 0 to 5, respectively, whereas the number of exons and introns in bHLH genes ranged from 23 to 1 and 0 to 22, respectively ( Figures 2 , 4 ). The majority of the genes (45 WRKY and 35 bHLH) from both families have three exons and two introns. Genes belonging to the same group of phylogenetic trees almost had parallel structures apart from a few genes. Among the SlWRKYs and SlbHLHs, Solyc07g005650 and Solyc09g063010 possess the longest structures, respectively. Among all SlbHLHs, a few genes have a complex structure, such as Solyc05g009220, Solyc10g079660, and Solyc03g114720. The exon–intron arrangement of the SlbHLH genes showed that there are exons lost or gained during the evolution process of the S. lycopersicum genome.

The MEME software was used to recognize conserved motifs to better understand the divergence and possible function of the WRKY and bHLH gene families ( Figures S1 , S2 ). Like gene structure, the motif distributions were similar within the same phylogenetic group. For instance, most of the members of group IV of the individual phylogenetic tree of SlWRKY have eight conserved motifs, and most of the members of group V of the SlbHLH individual tree have seven conserved motifs. The aqua blue-colored motifs were uniformly found in almost all of the WRKY and bHLH TFs, so these motifs, i.e., DPSIVITTYEGEHNH and KTDKASMLDEAINYIKELQKQ, significantly represented the conserved WRKY and bHLH domains, respectively. We did not identify a potentially conserved motif (aqua blue) in only two members (Solyc10g008250 and Solyc10g008270) of the bHLH family possibly due to lack of homology, sequence repetition, rearrangements, or disruption of alignment (Gribskov, 2019). These results of the gene structures and conserved motifs suggested that TFs clustered in the same group might be meaningful for gene evolution and their physiological roles.

We map the chromosomal locations of the SlWRKY and SlbHLH genes and illustrate that majority of the SlWRKYs and SlbHLHs are clustered at the chromosomal ends. The 81 SlWRKY genes are distributed in all S. lycopersicum chromosomes except chromosome 11. As represented in Figure 5 , most of the WRKY genes (19.7%) were located on chromosome 5, followed by chromosome 2 with 11.1% and chromosomes 3 and 8 with 9.8% of WRKY genes. In S. lycopersicum, out of 161, 160 bHLH genes were mapped from chromosomes 1 to 12, while the precise location of only one SlbHLH gene could not be determined. Chromosome 1 contains the highest number (13.7%) of bHLH genes, followed by chromosome 2 with 11.2% and chromosomes 3, 6, and 10 with 10% of bHLH genes, while chromosome 11 contains only 1.8% of the genes.

We further determined the tandem duplications of both WRKY and bHLH genes along the 12 S. lycopersicum chromosomes. Tandemly duplicated genes were characterized as an array of two or more homologous genes separated by 100 kb (Khan et al., 2021a). The gene duplication study showed that more than 50% of SlWRKY and SlbHLH genes were tandemly duplicated. As shown in Figure 6 , approximately 25 SlWRKY genes were tandemly duplicated unevenly on distinct chromosomes. Chromosome 2 had the most with eight SlWRKY gene pairs, followed by chromosome 5 with six SlWRKY gene pairs. Only one SlWRKY gene pair was discovered on chromosome 9, and chromosomes 8 and 11 had no WRKY gene pairs, while the Solyc05g045800 (SlWRKY67) and Solyc05g045710 (SlWRKY65) genes mapped on chromosome 5 had two pairs. The findings showed that segmental duplications assisted in the expression of SlWRKY genes. Similarly, approximately 28 SlbHLH gene pairs were discovered on distinct chromosomes in the S. lycopersicum genome. The maximum number of tandemly duplicated SlbHLH genes on chromosome 1 was 13, showing that a large number of SlWRKY genes on chromosome 1 were partly attributable to tandem gene duplication events.

Some WRKY, as well as bHLH proteins of S. lycopersicum, regulate the transcription process of the target genes when they form homodimer or heterodimer complexes to determine the DNA binding sites, and as a result, protein–protein interactions are critical in gene expression. The protein sequences in FASTA format were submitted to the STRING server. The protein–protein interaction network of differentially expressed proteins was constructed with the default setting. The retrieve included networks of the TFs (WRKY and bHLH), which highlight several hub proteins ( Table S2 ). WRKY proteins that have co-expression supported by a high confidence score (0.700) include WRKY70 (Solyc03g095770), which form a network with SlWRKY33A (Solyc06g066370), Solyc09g014990, Solyc03g116890, Solyc05g050300, and SlWRKY40 (Solyc06g068460), while Solyc05g014040 has co-expression with Solyc07g005650 and Solyc09g010960 ( Figure 7A ). In the bHLH family, the TFs with a higher co-expression score include Solyc07g043580–Solyc01g102300, Solyc08g075090–Solyc05g009220, and ICE1a (Solyc06g068870)–Solyc08g005050, while Solyc08g076820 has co-expression with Solyc05g053660 and Solyc09g091760 ( Figure 7B ).

To characterize the functions of SlWRKY and SlbHLH genes, GO enrichment analysis was conducted using the GO online tool. The results showed that the majority of the genes of both families were primarily involved in molecular functions and biological processes such as DNA binding (GO:0003677), regulation of DNA-templated transcription (GO:0006355), regulation of gene expression (GO:0010468), protein dimerization process (GO:0046983), regulation of RNA biosynthesis (GO:2001141), and nucleobase-containing compound biosynthetic process (GO:0034654) ( Figures 7C, D ).

WRKY and bHLH transcription factors play an important role in plants’ heavy metal stress defense (Guo et al., 2022). However, little information is available about the Cd-responsive WRKYs and bHLHs in tomato (S. lycopersicum). To determine the expression levels of the target genes in tomato plants under heavy metal stress, RNA-seq analysis was taken into consideration. The overall analysis of downregulated and upregulated genes was based on the differentially expressed genes in both Cd concentrations (1 and 2 mM of Cd) ( Figures 7E, F ). Genes were downregulated and upregulated based on log₂ fold-change >2 and log₂ fold-change <−2, respectively, with an adjusted p-value <0.05. Overall, 54 and 47 SlWRKY genes were differentially expressed in 1 and 2 mM of Cd stress, respectively, in which 43 and 40 genes were upregulated while 11 and 7 SlWRKY genes were found to be downregulated ( Figure 8A ). Similarly, approximately 46 common SlWRKYs were differentially expressed in both 1 and 2 mM of Cd stress ( Figure 8B ). On the other hand, a total of 70 and 69 SlbHLH genes were differentially expressed in 1 and 2 mM of Cd stress, respectively. However, most of the SlbHLH genes were to be downregulated in both 1 and 2 mM of Cd stress ( Figures 8A–C ). These findings suggest significant variations in both WRKY and bHLH genes between the control and Cd stress plants ( Figures 8D, E ). As shown in Figure 8D , our results revealed that SlWRKYs have almost the same expression levels under 1 and 2 mM of cadmium stress. Cadmium stress highly enhanced the expression of the following WRKY genes: Solyc05g007110 (SlWRKY76), Solyc02g094270 (SlWRKY38), Solyc08g067340 (SlWRKY46), Solyc06g048870 (SlWRKY19), Solyc09g014990 (SlWRKY33), Solyc02g021680 (SlWRKY35), Solyc08g067360 (SlWRKY45), Solyc04g051690 (SlWRKY51), and Solyc04g072070 (SlWRKY55). Contrastingly, cadmium stress downregulated the expression of Solyc01g104550 (SlWRKY9), Solyc10g084380 (SlWRKY44), Solyc09g010960 (SlWRKY49), Solyc08g081610 (SlWRKY29), Solyc04g051540 (SlWRKY13), and Solyc06g070990 (SlWRKY74). The results demonstrated that most of the SlbHLH genes’ overall expression level was downregulated during 1 and 2 mM of Cd stress, respectively ( Figure 8E ). The expression pattern of a few genes, including Solyc05g005300 (SlbHLH084), Solyc07g062200 (SlbHLH088), Solyc08g083170 (SlbHLH056), Solyc07g052930 (SlbHLH141), Solyc07g039570 (SlbHLH051), Solyc10g008260 (SlbHLH093), Solyc03g114230 (SlbHLH082), Solyc02g076920 (SlbHLH013), and Solyc01g106460 (SlbHLH007), was found to be upregulated. On the other hand, the expression of Solyc04g078690 (SlbHLH035), Solyc04g076240 (SlbHLH033), Solyc03g118310 (SlbHLH083), and Solyc08g076820 (SlbHLH146) was found to be significantly downregulated when S. lycopersicum plants were treated with cadmium which shows the importance of bHLH transcriptional factors in Cd stress tolerance.

To check the RNA-seq data accuracy, the relative expression levels of five SlWRKY (Solyc05g007110, Solyc02g094270, Solyc04g051690, Solyc01g104550, and Solyc10g084380) and five SlbHLH (Solyc07g062200, Solyc08g083170, Solyc01g106460, Solyc04g078690, and Solyc04g076240) genes were determined using qRT-PCR with specific primers ( Table S3 ). The expression patterns followed the same pattern as the RNA-seq data, and a substantial relationship was established between the RNA-seq and RT-qPCR results, demonstrating the reliability of RNA-seq data ( Figures 8F, G ).