The current research was conducted on the Tomato plant (S. lycopersicum c v. Yegwang). Seeds were surface sterilized according to previously reported protocols (Khan et al., 2023b). The seeds were properly washed with distilled water (dH₂O) thrice and then germinated at room temperature (28 ± 2°C; humidity 60 ± 5%) in total darkness for 2-3 days with the help of sterile filter paper. Seedlings were transferred to pots containing autoclaved soil which was composed of coco peat moss (10–15%), coco peat (45–50%), pertile (35–40), NH⁺ (ca.0.09 mg/g), KO (c.0.1 mg/g), zeolite (6-8%) with NO₃ (ca. 0.205 mg/g), and PO (ca.0.35 mg/g) in a greenhouse at a constant temperature (28 ± 2°C) under a 16/8 hr light and dark period, respectively. Three sets with 5 replicates per set were used for the experiment. Set 1 was kept as control, set 2 contained plants subjected to 1mM Cd stress, and set 3 contained plants subjected to 2mM Cd stress. The concentrations for Cd stress were selected based on the previous report (Khan et al., 2023a) and Cd stress was induced three days after transferring the plantlets into pots by soil irrigation method every three days. For the analyses, leaf samples from each replicate of all the sets were collected, and with the help of liquid nitrogen, the samples were immediately dried and stored at -80°C for RNA extraction.

Leaf samples freeze dried with liquid nitrogen were utilized to extract RNA according to previously reported protocols (Khan et al., 2023a). Briefly, The RNA obtained was diluted until the RNA quantity for each sample reached 100 ng. Three replicates from every treatment were carefully selected and sequenced for subsequent analysis. For RNA-seq analysis, the sequencing platform (Illumina HiSeq2000) was used as per the manufacturer’s procedure, resulting in 51-bp single-end reads. An efficient computational pipeline was used to find differences in gene regulation between the control and plants subjected to Cd stress. The quality was assessed through FastQC. Data was trimmed with the help of TrimGalore. To align and index the data to reference genome, HISAT2 was exploited. Feature Count (subread_v2.0.2) carried out read count quantification, and differential gene analysis was carried out with the help of DESeq2 in the R program. For DEG analysis, p-values were adjusted for multiple testing using the Benjamini-Hochberg procedure to control the false discovery rate (FDR). Genes were considered differentially expressed if they exhibited a log2 fold change (log2FC) greater than ±2, corresponding to a two-fold change in expression.

The S. lycopersicum and A. thaliana genes that encodes ATP binding cassettes (ABC) transporters were retrieved from Solgenomics (https://solgenomics.net accessed on 15 July 2023) and TAIR databases (https://www.arabidopsis.org/index.jsp accessed on 15 July 2023). The corresponding genome sequences and protein sequences were downloaded from Solanaceae Genomics Network (SGN, https://solgenomics.net/, accessed on 15 July 2023), and Arabidopsis Information Resource (TAIR, https://www.arabidopsis.org/index.jsp, accessed on 15 July 2023). Consequently, 154 ABC genes from S. lycopersicum and 130 ABC genes from Arabidopsis were identified. For valid identification and confirming the presence of the respective domains of the obtained protein sequences, the Simple Modular Architecture Research Tool (SMART, http://smart.embl.de/, accessed on 15 July 2023) was used at E-value <10-5. Only sequences containing the ABC transporter domain were retained, while sequences lacking this domain were excluded from further analysis. This validated dataset was then used for subsequent analysis

In order to classify the retrieved genes in S. lycopersicum (SlABCs), a phylogenetic tree was constructed by employing the Maximum likelihood (ML) technique with the help of software (MEGA 11), keeping the bootstrap replicates (BS) at 1000 times. The procedure above was repeated to construct a combined phylogenetic tree of 154 SlABCs with 130 ABCs from A. thaliana (AtABCs) as well, so as to determine the evolutionary relationship among the genes. The annotations of the resulting phylogenetic trees were performed using Version 1.4.4 of the FigTree (https://mybiosoftware.com/figtree-1-3-1-produce-figures-phylogenetic-trees.html#google_vignette).

Gene Structure Display Server (GSDS), an online tool (https://gsds.gao-lab.org/), was employed to investigate the gene structure, i.e., the number, structure, and organization of exons and introns in the identified SlABCs. In order to do so, the coding sequence (CDS) of SlABCs was aligned with the genomic sequence. For delving into the details of conserved Motifs in SlABCs, another online tool i.e. MEME (Multiple EM for Motif Elicitation) version 5.4.1 (https://meme-suite.org/meme/doc/meme.html) was employed by adjusting the motif number (maximum) to 20, keeping the optimum motif width between 6 to 200 amino acid (aa) residue and fixing the maximum repetitions to “any.”

The localization and distribution of SlABCs on the chromosomes of S. lycopersicum was accomplished by genome-wide collinearity analysis using TBtools software. The gene data required for this purpose was downloaded from the NCBI database.

The selected duplicate gene pairs were subjected to TBtools to compute the Ks and Ka values. The determined values were then used to determine the Ka/Ks ratio. The formula t = Ks/2r was used to calculate the divergence time with r being the representative value (1.5 × 10^-8) as reported in the literature (Ali et al., 2024).

An online tool called OrthoVenn2 (https://orthovenn2.bioinfotoolkits.net/home) was utilized for identification of orthologous groups among the ABC transporter encoding genes of S. lycopersicum and A. thaliana by analyzing their full-length sequences as reported by (Abiraami et al., 2023). For the mapping and analysis of Cis-regulatory modules, PlantCARE, an online website (http://bioinformatics.psb.ugent.be/webtools/plantcare/html) was explored using 2000 bp sequence upstream of SlABCs from S. lycopersicum genome with the help of a software (TBtools).

The protein sequences of 154 SlABCs were acquired in FASTA format and uploaded online to version 11.5 of STRING’s server (https://string-db.org/) to predict the nature of interaction among different proteins encoded by SlABCs and determine the putative function of respective proteins. The network construction was done with the help of selected genes (having significant confidence scores), whereas non-interactive genes were skimmed from the network.

For the prediction of annotated cellular components, molecular functioning and biological processes involving SlABCs, the enrichment analysis of Gene Ontology (GO) categories was performed with the help of online software agriGo (http://bioinfo.cau.edu.cn/agriGO/). Algorithm TopGo ‘elim’ was exploited for the said purpose using the loci of SlABCs, which were mined online from the SolGenomics website. Significantly enriched GO Terms (having False Discovery Rate (FDR) and P-value less than 0.05) were considered for the current study, as reported in literature (Khan et al., 2023b).

In order to authenticate the results of the RNA-Seq dataset, 10 genes were selected based on their significant differential expression observed under cadmium (Cd) stress conditions in the RNA-Seq analysis. These genes were chosen to represent the most biologically relevant responses to Cd stress. The selected genes were subjected to qRT-PCR. RNA dilution (0.1µg/µl final concentration) was carried out, followed by the cDNA production with the help of PCR BIOSYSTEMS’s qPCRBIO cDNA Synthesis Kits. qRT-PCR was conducted using BioFACTTM 2× Real-Time PCR Master Mix (Including SYBR^® Green I). A total of 20 μl sample containing master mix (10 μl), RNase free water (8 μl), cDNA, and primer (1 μl each) was utilized to carry out PCR reaction with the help of PCR machine (Step One Plus Real-Time PCR System, Life Technologies Holdings Pte Ltd., Singapore). The primers for selected genes were designed with the help of an online tool, i.e., Primer3 (https://bioinfo.ut.ee/primer3-0.4.0/), the details of which are available in ( Supplementary Table S1 ). A housekeeping gene, i.e., Actin, was used as an internal control for gene expression normalization. To analyze expression levels of selected genes from control plants and those subjected to Cd stress. Two technical and three biological repeats were employed for each sample in order to minimize experimental error. The conditions for the PCR reaction were as follows: 10 min at 94°C, followed by 35 cycles at 94°C (45 s), 65°C (45 s), and 72°C (1 min), with an extension step at 72°C (10 min). The gene amplification threshold was set at 0.1. Each sample was run three times with three different replicates.

The phylogenetic tree was constructed using the Maximum Likelihood (ML) method for in-depth analysis of the relationship between the SlABCs. As evident in Figure 1 , the results depicted that the identified SlABCs belong to 15 different groups, each labeled and designated by a specified color. More than half of the genes (55.8%) were grouped into 3 out of 15 groups, with Groups 5 and 6 having the maximum number of SlABCs (29 each), closely followed by Group 8 (28 SlABCs). The rest of the 12 groups collectively had 68 SlABCs (44.1%), with Group 11 having 13 SlABCs, Group 10 having 11 SlABCs, Group 1 having 9 SlABCs, Group 3 with 8 SlABCs, Group 13 with 5 SlABCs, Group 7, Group 9 and Group 14 having 4 SlABCs apiece and Group 2 and Group 15 with 3 SlABCs each. Groups 2 and 12 had the least SlABCs, i.e., 2 SlABCs each. A bootstrap value of over 99% in the majority of the clades validates the results above.

To further understand the evolutionary relationship of the SlABCs with AtABCs, a combined phylogenetic tree of 154 SlABCs and 130 ABC genes belonging to A. thaliana (AtABCs) was constructed using the aforementioned method. The resulting combined phylogenetic tree, as evident in Figure 2 , revealed 13 groups named A-M respectively. Apart from Group C and D, which only contained AtABCs (7 and 2 respectively), and L and M, which only contained SlABCs (5 and 2 respectively), the rest of the 9 groups contained ABC genes from both S. lycopersicum and A. thaliana. Group F contained the highest number of genes i.e. 102 (46 AtABCs and 56 SlABCs) followed by Group E with 57 (30 AtABCs and 27 SlABCs), Group G with 50 (18 AtABCs and 32 SlABCs), Group A with 23 (12 AtABCs and 11 SlABCs), Group I with 19 (7 AtABCs and 12 SlABCs), Group K with 9 (4 AtABCs and 5 SlABCs), Group B with 3 (2 AtABCs and 1 SlABC), Group J with 3 (1 AtABC and 2 SlABCs) and Group H with 2 (1 AtABC and 1 SlABC). Bootstrap value of over 95%-99% in all the clades validates the abovementioned results.

A total of 154 ABC genes were identified in S. lycopersicum, which were unequally distributed on chromosomes of S. lycopersicum ( Supplementary Table S2 ). Maximum number of SlABCs were present on chromosome 12 (20 SlABCs) followed by chromosome 6 (19 SlABCs), chromosome 3 and 11 (17 SlABCs each), chromosome 5 (15 SlABCs), chromosome 9 (12 SlABCs), chromosome 4 (11 SlABCs), chromosome 7 (10 SlABCs each), chromosome 2 and 8 (9 SlABCs each), chromosome 1 (7 SlABCs) and chromosome 10 (4 SlABCs). The precise location of 4 genes, however, was not identified. Apart from chromosomal distribution, other physical properties of the proteins encoded by SlABCs, like isoelectric point, molecular weight, alternative isoforms, number of exons and motifs, etc., also showed diversity ( Supplementary Table S2 ). The average pI value was recorded to be 7.64, with Solyc01g105400.2 showing the lowest pI value (4.31) and Solyc00g304030.1 displaying the highest value (10.43). Much variation was also noted in encoded proteins and gene length. The number of amino acids ranged from 53 to 1,909 whereas the gene length varied from 244 to 35,740 bp in 154 identified SlABCs. The predicted molecular weight ranged between 5.9 kDa and 213.17 kDa; the average molecular weight was 96.85 kDa. Subcellular localization revealed that ABC proteins were located in almost all major organelles, including plastid, vacuole, nucleus mitochondria, endoplasmic reticulum, chloroplast, and Golgi complex ( Supplementary Table S2 ).

A collinearity analysis was carried out to analyze the evolutionary relationship and duplication events in SlABCs. Duplication of genes is critical in increasing a gene family’s size. A number of SlABCs were found to be highly conserved on all the 12 chromosomes of S. lycopersicum suggesting multiple duplication events within the S. lycopersicum genome. The duplication included both intra-chromosomal and inter-chromosomal mutations ( Figure 3A ). Synteny analysis ( Figure 3B ) for identification of locus relationship among S. lycopersicum and A. thaliana genes showed the presence of 8 homologous genes. These genes were located on 7 out of the 12 chromosomes of S. lycopersicum. Interestingly, the gene on chromosome 6 of S. lycopersicum had homologous genes on chromosome 1. 2. 3, and 4 of A. thaliana, whereas chromosome 12 had 2 different genes having homologous genes on chromosome 3 and 4 of A. thaliana.

In order to estimate the balance between different types of mutations, the Ks and Ka values of selected homologous gene pairs were calculated. The values were used to calculate the ratio between synonymous and non-synonymous values, the mean for which was 0.411. The results showed that the homologous genes evolved through purifying selection, suggesting that mutations that could deteriorate the genes were removed. The whole process was estimated to have occurred approximately 9 to 79 mya depending upon the homologous gene pair in question ( Table 1 ).

Examining gene structure and motif composition offers valuable insights into evolutionary history and functional diversity. To explore these aspects, the conserved motifs and exon-intron organization of the 154 SlABC genes in S. lycopersicum were analyzed ( Figures 4 , 5 ). The number of exons in SlABC genes ranged from 1 to 40, while introns varied between 0 and 39, indicating considerable structural diversity within this gene family. MEME software identified 20 distinct conserved motifs, reflecting potential functional specialization among SlABC genes. Among these motifs, several were identified as Nucleotide Binding Domains (NBDs), crucial for ATP binding and hydrolysis. Specifically, Motif 2 contained the Walker A motif, a key ATP-binding site, while Motif 5 harbored an ABC transporter signature motif essential for transporter function. Additionally, Motif 12 was associated with the Walker B motif, which plays a fundamental role in ATP hydrolysis. Motif 13 was identified as another ATP-binding region commonly found in ABC transporters. The analysis further revealed that full-length ABC transporters predominantly contained two NBDs, whereas half-size ABC transporters had only one. The position of these domains varied across different genes, appearing at either the N-terminal, C-terminal, or in a duplicated arrangement, depending on the specific subfamily.

The conservation of these motifs across multiple SlABC genes highlights their strong functional preservation, particularly in ATP-dependent substrate transport. These findings suggest that despite sequence variations, the core functional components of the SlABC gene family remain highly conserved.

A total of 42 ortholog clusters (1-42) were identified which are given in Supplementary Figure S1 . The orange dots represent A. thaliana ABC genes (AtABCs), whereas the blue dots represent SlABCs. The 42 designated clusters show different configurations, including many to many (eight clusters i.e. 1-7 and 10), many to one (six clusters i.e., 8-9, 13-16), one to many (5 clusters i.e. 11-12, 17-19) and one to one (23 clusters, i.e. 20-42) specific ABC members. Interestingly, the orthologs belonging to the same cluster had the same number of motifs, For example, all genes in cluster 2 had the same number of motifs (16) and a similar motif sequence. Similarly, genes belonging to clusters 11 and 12 had seven motifs with a similar motif sequence. Similarly, the number of exons and introns within the same cluster was also similar. For example, all the genes in Cluster 4 had 18 exons and 17 introns, whereas genes in Cluster 5 had a single exon with no introns. The number and sequence of SlABC genes in clusters 4 and 5, respectively, are also similar. The same pattern was observed for most of the remaining clusters as well. The above discussion validates the fact that the orthologs have evolved from a common ancestral gene.

According to Figure 2 , the combined phylogenetic tree divided SlABCs and AtABCs into 13 groups. Except for Group B, which did not have any orthologs, all the other groups have a number of orthologs i.e. Group F contains 19 orthologs, Group E contained a total of 13 orthologs, Group G have six orthologs, Group A and Group I have five orthologs each, Group K contains two orthologs and Group J and Group H have one ortholog each. Although A. thaliana and S. lycopersicum belong to different families that diverged millions of years ago, most of the ABCs have orthologous relationships, suggesting that the nature and functions of AtABCs and SlABCs would show a similar pattern.

Various hormones play a critical role in enabling the plant to combat different biotic and abiotic stress conditions. In this regard, 137 SlABCs were found to play a role in hormone-responsive cis-regulatory elements ( Figure 6 ). A total of 1150 cis-regulatory modules comprising 15 different elements were identified in which ABRE, an ABA-responsive element, was the most repeated cis-regulatory element present (279) found in 114 SlABCs followed by CGTCA-motif (155), a Methyl jasmonate (MeJA) responsive element and TGACG-motif (155) which is involved in the regulation of salicylic acid (SA) and MeJA signaling pathways which were identified in 54 and 53 SlABCs respectively. The number of cis-regulatory modules found in individual SlABCs showed tremendous diversity and ranged between 1-24 with Solyc12g098870.1 having the most number of cis-regulatory modules (7 TGACG-motifs and CGTCA-motifs, 6 ABRE and 1 O₂-site, TC-rich repeats and TATC-box) and two SlABCs (Solyc11g018680.1 and Solyc09g042280.1) having only one cis-regulatory module. Interestingly, only one SlABC (Solyc11g065350.1) contained the SARE motif, a SA-responsive element. From the above discussion, it is evident that all of the SlABCs contain at least one stress-responsive cis-regulatory element, suggesting their role in combating a variety of biotic and abiotic stress to which plant is constantly exposed due to their sessile nature.

Protein-protein interactions (PPI) offer valuable insights into the biological functions of proteins, as those involved in similar pathways often interact with each other. To better understand the functional associations of SlABC transporters, a protein-protein interaction (PPI) network was constructed using STRING ( Figure 7A ). The analysis revealed a well-connected network, where several SlABC transporters formed central hubs, suggesting their significant role in intracellular transport and regulatory processes. The functional enrichment of interacting proteins highlighted key biological processes, including ATP-dependent transport, xenobiotic efflux, lipid transport, and long-chain fatty acid metabolism. Notably, the presence of ABC-type xenobiotic transporters suggests that some SlABC proteins may be involved in detoxification, while interactions with lipid transporters indicate a potential role in maintaining membrane stability and homeostasis. Gene ontology (GO) enrichment further confirmed associations with transmembrane transport, ATP binding, and efflux mechanisms, emphasizing their importance in energy-dependent substrate transport.

These findings prove that SlABC transporters are functionally diverse and play crucial roles in cellular transport, stress adaptation, and metabolic processes in S. lycopersicum. The PPI network suggests that certain SlABC transporters act as key regulators of nutrient transport and environmental stress responses, potentially contributing to overall plant resilience and development.

To predict the annotated function of the identified SlABCs, Gene Ontology (GO) enrichment analysis was carried out through online GO tools. As evident in Figure 7 , the SlABCs were involved in a number of biological ( Figure 7B ), molecular ( Figure 7C ), and cellular components ( Figure 7D ). Key biological processes include transport (GO:0006810), localization (GO:0051179), the establishment of localization (GO:0051234), and transmembrane transport (GO:0055085), whereas molecular functions included a wide range of activities like active transmembrane transporter activity (GO:0022804), ATP binding (GO:0005524), adenyl ribonucleotide binding (GO:0032559), purine ribonucleoside triphosphate binding (GO:0035639), transmembrane transporter activity (GO:0022857), ribonucleotide binding (GO:0032553), carbohydrate derivative binding (GO:0097367), anion binding (GO:0043168) and ABC-type transporter activity (GO:0140359) etc.

Plants in their natural habitat had to face a number of extreme conditions like drought, salinity, and heavy metal-induced stress. A variety of ABC isoforms are well known for their role in the reclamation of heavy metal-induced stress; however, the response of SlABCs to Cd-induced stress has not been explored to the best of our knowledge. In order to address the existing research gap, the S. lycopersicum plantlets were exposed to two different concentrations of Cd (1 mM and 2 mM) and expression levels of targeted SlABCs was evaluated using RNA-seq analysis ( Figure 8 ). About 105 selected SlABCs were found to be expressed during Cd stress, with several genes showing differential expression patterns. 15 SlABCs were expressed in response to 1mM Cd stress alone, 8 SlABCs were expressed in response to 2 mM Cd, and 50 SlABCs showed expression in both 1mM and 2mM Cd stress ( Figure 8B ). In the case of 1mM Cd stress, 10 SlABCs were upregulated, whereas 5 SlABCs were downregulated. Similarly, 3 SlABCs were upregulated, and 5 SlABCs were downregulated in response to 2mM Cd stress, whereas 35 SlABCs were upregulated and 15 SlABCs were downregulated in response to both 1mM and 2mM Cd stress ( Figure 8C ).

Using specific primers, the RNA-seq-based gene expression data was validated by subjecting 10 SlABCs, i.e., Solyc03g007690.1, Solyc04g010200.1, Solyc04g070970.4, Solyc05g053610.2, Solyc05g056470.1, Solyc06g036240.2, Solyc06g075020.3, Solyc07g063400.3, Solyc08g067610.3, and Solyc08g067620.2 to qRT-PCR for analyzing the relative gene expression levels. There was a strong correlation between gene expression results obtained by both RT-PCR and qRT-PCR. For instance, 5 SlABCs i.e. Solyc03g007690.1, Solyc04g010200.1, Solyc04g070970.4, Solyc05g056470.1 and Solyc06g075020.3 were downregulated in both RNA-seq ( Figure 8D ) and qRT-PCR ( Figure 9 ) data whereas the other 5 SlABCs i.e. Solyc05g053610.2, Solyc06g036240.2, Solyc07g063400.3, Solyc08g067610.3, and Solyc08g067620.2 were all upregulated in both RNA-seq and qRT-PCR based gene expression analysis thereby validating the accuracy our results.