Whole genome of tomato (S. lycopercicum) was downloaded from Solanaceae Genomics Network (https://solgenomics.net/) database to construct a local database by employing BioEide7.0 (bioedit.software.informer.com). In order to retrieve tomato SlHsp70 gene sequences, Arabidopsis thaliana AtHsp70 gene sequences were BLAST-analyzed against newly constrctued local tomato database. The Hidden Markov Model (HMM) was employed to validate SlHsp70 gene sequences and their functional domains particularly by running against each AtHsp70 gene (Pfam: PF00012). Incomplete raw reads were excluded and assigned specific name to each gene according to Mendel database for plant gene families listed in Commission on Plant Gene Nomenclature (CPGN) (mbclserver.rutgers.edu/CPGN/), International Society of Plant Molecular Biology (ISPMB) [32]. In order to further validate retrieved SlHsp70 genes, BLAST analyses were also performed in NCBI genome database by selecting tomato whole genome (blast.ncbi.nlm.nih.gov/blast.cgi), SPud DB tomato Solanaceae Genomics Network (solgenomics.net) and phytozome (phytozome.jgi.doe.gov/).

All retrieved SlHSP70 polypeptide sequences were further validated at E-value 10${}^{-5}$ for detection of functional domains by deploying SMART (smart.embl-heidelberg.de/) tools. All SlHsp70 gene sequences were analysed for chromosomal localization. SlHSP70 family proteins were in-silico characterized for their molecular weights (MW), total number of atoms, total number of amino acids, Aliphatic index, isoelectric point (pI) and instability index by employing EXPASY PROTOPARAM (expasy.org/tools/protparam.html) [33, 34]. Theoretical pI and total MW of each polypeptide was analyzed by employing ProtParam (web.expasy.org/portparam).

Arabidopsis AtHSP70 family proteins were downloaded from Arabidopsis Information Resource (TAIR) (arabidopsis.org) and potato StHSP70 family proteins from the Spud DB Potato Genomics Resources (solanaceae.plantbiology.msu.edu/). Multiple sequence allignments of HSP70 polypeptide sequences of Arabidopsis, potato and tomato were performed in ClustalX 2.0 with default parameters [35]. Phylogenetic analysis were performed to construct an outgroup rooted tree by using 20 A. thaliana AtHSP70, 19 S. tuberosum StHSP70 and 23 S. lycopercicum SlHSP70 protein sequences by performing alignment using Neighbor-Joining (NJ) method in MEGA7.0 with following parameters: 2000 bootstrap replication values, pair-wise gap deletion mode and Poisson model [36, 37].

The phytozome plant genome database (phytozome.jgi.doe.gov/) was explored by selecting tomato to collect information obout localization of SlHsp70 genes and a genetic map was constructed by employing MapChart tool [38]. Genes of similar species were placed in the same group of an outgroup rooted tree, defined as co-paralogs, which were further analyzed to identify tandem duplicatons, segmental duplications and their coordinates. To investigate evolutionary events such as gene duplicaiton and divergence for vigourous gene functions as well as genetic expansion, we explored PGDD (Plant Genome Duplication Database) by deploying circos. Co-paralogs were considered tandemly duplicated when both of them were non-homologous and distance between their loci was 100 kb [39] and segmentally duplicated if located on different chromosomes [40]. Smith-Waterman algorithm (ebi.ac.uk/Tools/psa/) was employed to calculate local alignment of two proteins and synteny analyses were performed by employing circos with default parameters (circos.ca/). Finally, amino acid sequences of all SlHSP70s were used to construct protein-protein interaction network with the help of STRING (version: 11.5) proteome database (string-db.org/) by selecting Solanum lycopercicum as an organism.

Genomic and CDS sequences of each SlHsp70 gene were analyzed to examine number and order of introns and exons with the help of Genes Structure Display Server 2.0 (GSDS 2.0) program (gsds.cbi.pku.edu.cn/index.php) by adjusting following parameters: (a) Intron: color-black, shape-line and line width-3, (b) CDS: color-dark maroon, shape-round corner rectangle and height-12 and (c) UTR: color-blue, shape-rectangle and height-10. Protein sequence of each SlHSP70 was also analyzed to identify conserved motifs with the help of Multiple EM for motif elicitation (MEME) tool (meme.nbcr.net/meme3/meme.html) with following parameters: number of motifs-10 and optimum amino acid residues per motif were 6 to 200 [41].

Phyre2 web (sbg.bio.ic.ac.uk/phyre2/) was utilised in intensive mode to undertake protein modelling of SlHSP70 polypeptides [42]. All predicted models of SlHSP70 proteins were based on c5tkyA, c5e84B, c3d2fC, c2khoA, c3c7nB and c5obuA templates with 100% identification. In-silico direct physical and indirect functional interactions among proteins were predicted by using STRING (version 11.5) database. GO enrichment analyses were performed to predict localization of SlHSP70 proteins within cellular components (red bars), expected molecular function (green bars) and participation in possbile biological processes (blue bars) (Fig. 1 & Supplementary Table 1). To perform GO enrichment analysis for each SlHSP70 protein, OmicsBox and Blast2GO v3.0.11 (www.blast2go.com) were employed [43].

Fig. 1.

GO enrichment analysis of SlHsp70 genes in S. lycopercicum. Red columns represent cellular components in which SlHsp70 genes displayed expression, blue columns represent biological processes in which SlHsp70 genes played their specific roles and green columns represent expected molecular function of SlHsp70 genes.

We analyzed RNA-seq data under normal conditions and constructed a heatmap to explore expression profiles of all SlHsp70 family genes in different tissues of S. lycopercium including flower bud, unopened flowers, fully opened flowers, fruits (1 cm, 2 cm and 3 cm), mature green fruits, breaker fruits, breaker + 10 fruits, leaves and roots (Fig. 2 and Supplementary Table 2). RNA-seq data and PFKM values were downloaded from tomato functional genomics database (ted.bti.cornell.edu/pgsc_download.shtml) and analyzed using cufflinks v2.2.1. FPKM values were divided by their mean, transformed into log2 ratio and clustered into expression data in the form of heat map (heatmapper.ca/) with the help of MeV4.5 with default parameters [44, 45].

Fig. 2.

Expression level of of all 23 SlHsp70 in different tissues including roots, leaves, buds, flowers and fruits of different stages of S. lycopercicum based on RNA-seq (http://ted.bti.cornell.edu/).

The seeds of tomato (line M82) were cultured in greenhouse of Yibin University in mid Autumn 2020. In order to perform surface sterilization, seeds were soaked in 10% hypochlorous acid for 5 mins and then washed thrice with ddH${}_2$O. After sterilization, seeds were placed on wet filter paper for germination, germinated seedlings were detached from paper and anchored in pots filled with Miracle-Gro® and Metro-Mix® 200 soil mixture and placed in greenhouse under following conditions 16 hours light at 27 ºC, 8 hours dark at 18 ºC and 70% humidity. Thirty days old tomato seedlings were uprooted, roots were washed with ddH${}_2$O to remove soil particles and dipped in 1/2 Hoagland solution (pH 6.0) supplemented with different concentrations of five trace metal elements such as 1 mM MnSO${}_4$, 0.1 mM CdCl${}_2$, 0.5 mM FeSO${}_4$, 0.1 mM CoCl${}_2$, 0.5 mM ZnSO${}_4$ and normal 1/2 Hoagland solution as control (CK) [46, 47]. To extract RNA for qRT-PCR anslysis, leaves and roots of tomato seedlings treated with Hoagland solution supplemented with different concentrations of trace metal elements and normal Hoagland solution was removed after 24 h and stored in liquid nitrogen.

We used TRIzol™ Reagent (Thermo Fisher Scientific, USA) for total RNA extraction and SuperMix Kit (TransGen, Beijing) was used for cDNA synthesis. Gene specific primers were manually designed (Supplementary Table 3) and $\beta$-actin was selected as an internal control during qRT-PCR [48]. qRT-PCR was performed with following recation mixture 10 $\mu$L of SYBR premix Taq 2X, 1 $\mu$L of cDNA as a template, 0.5 $\mu$L of each forward primer and reverse primer and 8 $\mu$L of ddH${}_2$O. PCR conditions were as follow: 95 °C for 10 min, 95 °C for 15 s, total 40 cycles, 60 °C for 60 s and three scientific replications. Relative expression level of each SlHsp70 gene was measured with the help of Livak equation (2${}^{-\mathrm{\Delta }\mathrm{\Delta }\text{CT}}$) [49].

Hsp70 genes sequences were retrieved by downloading and subsequently analyzing S. lycopercicum genome database, incomplete raw reads were excluded and finally 23 complete candidate Hsp70 gene sequences were selected. All genes were assigned scientific name given as SlHsp70-1 to SlHsp70-23 (Table 1). Gene synteny analysis revealed that all 12 tomoto chromosomes harbored Hsp70 genes except chromosome 5. Total number of amino acids in each SlHSP70 polypeptide were 186 to 890, and molecular weight (MW) was 21273.47 Da to 98787.40 Da (Table 1). Comparatively, high contents of negative and acidic amino acids such as L-glutamic acid and L-alpha-aspartyl residues (Asp + Glu) were observed in almost all SlHSP70 proteins except SlHSP70-12 (Table 1). Isoelectric point (pI) of all 23 genes was also acidic (Table 1).

Phylogenetic analysis revealed distribution of all SlHsp70 genes in four groups I-IV (Fig. 3). Highest number of Hsp70 genes were placed in Group I of phylogenetic tree which comprised of 6, 9 and 8 Hsp70 genes of Arabidopsis, potato and tomato, respectively. Group II was smallest which comprised of five tomato SlHsp70 genes, two potato Hsp70 genes and three Arabidopsis AtHsp70 genes, respecitvely. Group III comprised of four tomato SlHsp70 genes, three potato Hsp70 genes and five Arabidopsis AtHsp70 genes. Group IV was second largest group which compried of six tomato SlHsp70 genes, five potato Hsp70 genes and five Arabidopsis AtHsp70 genes (Fig. 3).

Fig. 3.

Phylogenetic tree comprised of 23 S. lycopercicum (red circle), 19 Arabidopsis thaliana (green circle) and 19 S. tuberosum (purple circle) HSP70 protein sequences. ClustalX 2.0 was employed for protein alignment. Neighbor-Joining (NJ) method was used at 2000 bootstrap value to construct a phylogenetic tree in MEGA 7.0.

An uneven distribution of all SlHsp70 genes on different chromosomes was detected except chromosome number 5 which did not harbor any SlHsp70 gene. The highest number of SlHsp70 family genes localized on any chromosome were following four: SlHsp70-1, SlHsp70-5, SlHsp70-10 and SlHsp70-21, which were present on chromosome 1 (Fig. 4). Both segment and tandem duplication revealed 120 collinearity gene pairs with 50–100% duplication events (Supplementary Table 4). We noted 4 sister pairs among all 23 SlHsp70 genes. Multiple segment duplication events of SlHsp70-3, SlHsp70-10 and SlHsp70-13 genes were also observed. On the contrary, four parlog pairs were observed having distance less then 5 kb which were tandemly duplicated gene clusters and localized on chromosome 1, 3, 8, 11, and 12 (Fig. 4 and Supplementary Table 5).

Fig. 4.

Gene synteny analysis of SlHsp70 genes in S. lycopercicum. Blue lines represent orthologs and paralogs to express segmental duplication while red boxes represent tandem duplication.

On the base of structure, all 23 SlHsp70 gene family members were sub-divided into A, B, C, D and E sub-families (Fig. 5a). Largest subfamily was A containing 8 genes, subfamily E was second largest with 6 genes and subfamily D was smallest with only 1 gene (Fig. 5a). Except for SlHsp70-9 and -13, which had no introns, the number of exons and introns varied across all sub-families. Highest degree of similarity in exon and intron numbers was observed with in sub-families (Fig. 5c). Each motif of SlHPS70 proteins comprised of 50 amino acids except motif 4, 6, and 8 which contained only 41, 33 and 29 amino acids, respectively. The motifs 1, 3 and 8 were found in all subfamilies except D, and were followed by motifs 4, 6 and 7, which were likewise found in all subfamilies but two members of E. Noticeably, number, order and types of motifs were similar within a subfamily but different among subfamilies (Fig. 5b and Table 2).

Fig. 5.

Phylogenetic analysis, gene structure analysis and conserved motifs of SlHsp70 genes in S. lycopercicum. (a) The neighbor-joining (NJ) method was used at 2000 bootstrap value to construct a phylogenetic tree of SlHSP70 amino acid sequences in MEGA 7.0. (b) Ten conserved motifs of all SlHSP70 proteins were presented in unique colour symbol shown in box. (c) Dark lines represent exons, black lines represent introns annd blue lines represent UTRs. Specific size scale is given below.

Table 2.Conserved motifs of 10 selected SlHSP70 proteins in S. lycopercicum.

Motif	Logo	Most probable matches	Width
1		VKBAVVTVPAYFNDSQRQATKDAGVIAGLNVLRIINEPTAAAJAYGLDKK	50
2		LLDVTPLSLGJETAGGVMTKLIPRNTTIPTKKEQVFSTYSDNQPGVLIQV	50
3		EKNVLVFDLGGGTFDVSJLTIEEGIFEVKATAGDTHLGGEDFDNRLVNHF	50
4		TRARFEELNMDLFRKCMEPVEKCLRDAKLDKSDIHEVVLVGGSTRIPKVQ	41
5		EGERARTKDNNLLGKFELSGIPPAPRGVPQIEVCFDIDANGILNVSAEDK	50
6		FNGKEPCKSINPDEAVAYGAAVQAAILSG	33
7		ERLIGDAAKNQAAMNPENTVFDAKRLIGRRFSDP	50
8		FKRKHKKDISGBPRALRRLRTACERAKRTLSSTAQTTIEIDSLYEGIDFY	29
9		FKRKHKKDISGBPRALRRLRTACERAKRTLSSTAQTTIEIDSLYEGIDFY	21
10		YKGEEKQFSPEEISAMVLTKMKEIAEAFL	29

Amino acid sequences of SlHSP70 were analyzed in-silico to predict 3D protein structures as 3D conformation guarantees specific function of any protein (Fig. 6 & Supplementary Table 6). Except for SlHSP70-1 and -5, the rest of the SlHSP70 proteins were modelled using the c3d2fC template at a confidence level of 100 percent., while modeling of SlHSP70-12, -13 and -20 proteins was performed using c2khoA, c5gjjA, c5tkyA and C5nnrD templates. STRING database analysis revealed 22 nodes, 97 edges and following 13 local network clusters: 4141, 4411, 4655, 4661, 4020, 4021, 4023, 4025, 4084, 3997, 4339, 4682 and 4138. Among all, the biggest cluster was 3997 which comprised of 12 SlHsp70 proteins (Supplementary Table 7). Protein analysis also revealed existence of two common protien domains (PF00012 in all 23 SlHSP70 proteins and PF06723 in 19 SlHSP70 proteins) and five KEGG pathways (such as sly04141 in 14 SlHSP70 proteins, sly03060 in 3 SlHSP70 proteins, sly04144 in 10 SlHSP70 proteins, sly03040 in 10 SlHSP70 proteins and sly03018 in 2 proteins) (Supplementary Table 8). Sub-cellular distribution percentages of SlHSP70 proteins were 2/4% in vaculoar membranes, 3.5/8% in chloroplast and 9/42% in endoplasmic reticulum (ER). SlHsp70-8 was localized in 21/28 sub-cellular compartments. Overall percentages of SlHSP70 proteins in different biological processes were as follow: heat and cadmium resistance response was 3/4%, cellular stress response was 17/23% and chemical stress response was 15/23%.

Fig. 6.

Prediction of 3D structures of SlHSP70 proteins in S. lycopercicum. Phyre 2 server was used in an intensive mode to generat protein models, visualized by rainbow colours in direction from N to C terminus.

Tissue specific expression of the gene SlHsp70-1, -5, -9, -10, -13, -14, -15, -16, -17, -18, -20 and -21 diplayed upregulated expression level in breaker + 10 fruits while SlHsp70-8 exhibited highest expression level in breaker fruits. SlHsp70-11 displayed upregulated expression level in fully opened flowers and mature green fruits while SlHsp70-22 showed mild expression level in 2 cm fruit and upregulated expression level in fully opened flowers. Only single gene with upregulated expression level in leaves was SlHsp70-12 but it also displayed very mild expression in unopened flower bud. SlHsp70-6 exhibited mild expression level in roots, 2 cm fruit, breaker + 10 fruits and mature green fruits. SlHsp70-19 displayed mild expression in roots, 1 cm fruit, 2 cm fruit and unopened flower bud. Finally, SlHsp70-23 showed mild expression in mature green fruit and 3 cm fruit (Fig. 2).

Each SlHsp70 gene exhibited differential expression level on treatment of different trace metal elements in both leaves and roots. The roots of plants treated with Cd${}^{2+}$ and Co${}^{2+}$ solutions displayed upregulated expression of all SlHsp70 genes except SlHsp70-7, -17 and -20, while SlHsp70-11 specifically showed upregulated expression in roots (Fig. 7). Similarly, Cd${}^{2+}$ and Co${}^{2+}$ solutions treatment resulted in upregulation of expression of all SlHsp70 genes in leaves except SlHsp70-7, -11, -17 and -20. Fe${}^{2+}$ treatment resulted in increased expression of SlHsp70-17 and -23 but decreased expression of SlHsp70-18, -20 and -21. Mn${}^{2+}$ treatment resulted in decreased expression of SlHsp70-14 and SlHsp70-21, but increased expression of SlHsp70-3, -8, -12 and -23. Zn${}^{2+}$ treatment resluted in decresed expression of SlHsp70-7, -14, -18 and -21 but increased expression of SlHsp70-1, -3, -8, -10, -12 and -23 (Fig. 7).

Fig. 7.

qRT-PCR analysis of all 23 SlHsp70 genes in root and leaf tissues under different metal stresses. $\beta$-actin gene was employed as an internal control. Standard deviation among three scientific repeats was shown by drawing error bars. Mean expression of three biological replicates of each SlHsp70 under Cd, Co, Fe, Mn and Zn stress was taken by deploying t-test at p 0.05 significance (n = 9). CK means control and asterisk (*) represents significant differential expression of SlHsp70 genes in root and leaf under stress of specific heavy metal.