Genome-wide identification of C2H2 zinc-finger genes and their expression patterns under heat stress in tomato (Solanum lycopersicum L.)

Genome-wide analysis of the C2H2-ZF family in tomato

The tomato C2H2-ZFP family members were identified using their homology to the Arabidopsis thaliana C2H2-ZFP sequences from the TAIR10 database (File S1). The tomato genome (The Tomato Genome, 2012) and protein sequences (https://solgenomics.net/organism/Solanum_lycopersicum/genome) were used to construct a local protein database using Geneious v4.8.5 software (http://www.geneious.com/; Biomatters, Auckland, New Zealand), which was then BLASTP searched using the sequences of the A. thaliana C2H2-ZFPs (obtained from the Arabidopsis Information Resource; TAIR10, https://www.arabidopsis.org) as queries, with an E-value cut-off ≤1 × 10⁻²⁰. Subsequently, the Hidden Markov Model (HMM) profiles of the C2H2-ZFP sequences (Pfam ID: PF00096) were downloaded from the Pfam database (http://pfam.xfam.org) and used to validate the identity of all candidate C2H2-ZF gene members.

All the obtained C2H2-ZFP sequences were further confirmed using the NCBI Conserved Domain Database (CDD; https://www.ncbi.nlm.nih.gov/cdd/) with the default parameters. Proteins without C2H2-ZF domains were removed. The locus ID and chromosomal location information of each tomato C2H2-ZF gene family member were obtained from the genome annotation file (File S1), and the lengths of the coding sequences (CDSs) were determined by performing BLASTn searches against the Sol Genomics Network database (https://solgenomics.net). The physicochemical properties of the deduced proteins, including the molecular weight (MW), isoelectric point (pI), and grand average of hydropathy (GRAVY) values, were determined using the ExPASy-ProtParam tool (http://web.expasy.org/protparam/).

Phylogenetic analysis and gene duplication

To identify the evolutionary relationships of the C2H2-ZF gene family members, all C2H2-ZFP sequences were aligned using ClustalW2 software under the default settings (Larkin et al., 2007). The aligned sequences were then subjected to a phylogenetic analysis using MEGA v6.0 (Tokyo Metropolitan University, Tokyo, Japan; Tamura et al., 2013). Subsequently, we compared the phylogenetic trees constructed and tested by different methods, including Maximum Likelihood (ML), Neighbor Joining (NJ), unweighted pair-group method with arithmetic means (UPGMA) and Maximum Parsimony (MP) methods, respectively. The phylogenetic trees were constructed using different methods with 1,000 replicate bootstrap tests. The ML trees was calculated using the ProtML program under the JTT model (Guo et al., 2008), NJ trees was obtained using the JTT+I+G substitution model (Pan et al., 2018), UPGMA and MP trees were generated in MEGA v6.0 (Tamura et al., 2013) with the default parameters. Finally, the phylogenetic trees were visualized using FigTree v1.4.2 (http://tree.bio.ed.ac.uk/software/figtree/).

Based on the GFF genome files from the Sol Genomics Network database (https://solgenomics.net; The Tomato Genome, 2012), a map of the distribution of the C2H2-ZF genes was drawn for each chromosome using MapChart v2.0 (Voorrips, 2002).

Gene structure, conserved motif and functional annotation analyses

Genome DNA and the corresponding CDSs of the putative C2H2-ZF genes were obtained from the Sol Genomics Network database (https://solgenomics.net; The Tomato Genome, 2012), then analyzed using the Gene Structure Display Server (GSDS v2.0; http://gsds.cbi.pku.edu.cn/index.php) to obtain information on the exon/intron structures. MEME v5.0.3 (http://meme-suite.org/tools/meme) was used to predict the corresponding conserved motifs (Bailey et al., 2009) with the following parameters: optimum motif widths of 6–300 residues, any repetition, and a maximum of 10 motifs (Pan et al., 2018). Each motif with an E-value <1 × 10⁻¹⁰ was retained for motif detection. In addition, gene ontology analysis of tomato C2H2 family genes was performed using the Blast2GO program (Conesa et al., 2005) with the default parameters.

Expression analysis of the C2H2-ZF family members

The expression profiles of the C2H2-ZF genes were measured using the publicly available tomato RNA-Seq datasets retrieved from the TomExpress database v17.0.0 (Available online: http://tomexpress.toulouse.inra.fr/; Zouine et al., 2017). Subsequently, the expression levels had been normalized in tomato as the published method (Maza et al., 2013), and the pipeline was described below. The reads were mapped to the tomato genomes SL2.40 by the TopHat 2.0.0 software with default parameters (The Tomato Genome, 2012; Trapnell et al., 2012). Gene expression levels were assessed with Cufflinks software with default parameters and using the tomato gene annotation file ITAG3.2 (The Tomato Genome, 2012; Trapnell et al., 2012). Eventually, the visualized heatmaps were generated using Heatmap Illustrator 1.0.

Plant materials and heat treatment

The wild-type tomato (Solanum lycopersicum cv. Ailsa Craig) was grown in 15-cm pots containing a soil mix of humus, vermiculite, and perlite in a ratio of 3:1:1. The plants were grown under standard glasshouse conditions (16 h light and 8 h dark at 25  ± 2 °C). After one month of cultivation, the plants were subjected to a heat stress condition (42 °C) for 12 h or 24 h, while those collected at 0 h were used as the controls. All samples (leaves, stems, and roots) were collected and immediately frozen in liquid nitrogen for RNA extraction. Three biological repetitions were performed for each treatment.

RNA isolation and qRT-PCR analysis

Total RNA was isolated from the roots, stems, and leaves using the RNAiso Plus Kit (TaKaRa Bio, Dalian, China). 1-µg aliquot of RNA was treated with DNase I (Takara Bio, Dalian, China) and transcribed into cDNA using the Oligo dT-Adaptor Primers and an RNA PCR Kit (AMV) v3.0 (Takara Bio, Dalian, China). The qRT-PCR analysis was performed on a CFX96 Real-Time PCR system (Bio-Rad Laboratories, Hercules, CA, USA) using Eva Green SMA (Bio-Rad Laboratories). To evaluate relative gene expression levels, the SlEF1-α gene (Solyc06g005060), SlACT (Solyc03g078400), and SlUBI3 (Solyc01g056940) were used as the internal reference, and 35 pairs of gene-specific primers were used for the qRT-PCR (Table S1). The relative expression levels of these genes were calculated using the 2^−ΔΔCT method (Pan et al., 2012). To ensure the statistical credibility of the results, all experiments were performed with three biological replicates and three technical replicates. And data were compared using t-test.

Genome-wide identification and characterization of the C2H2-ZFs in tomato

Using BLAST and the HMM profiles of the C2H2-ZF domain, a total of 116 tomato cDNAs in the tomato genome (cDNA release 3.20; https://solgenomics.net/organism/Solanum_lycopersicum/genome) were annotated as encoding C2H2-ZFPs. All the C2H2-ZF candidates were analyzed using the NCBI Conserved Domain Database to verify the presence of the C2H2-ZF domain. Finally, 104 C2H2-ZFPs were confirmed in tomato (File S1). Correspondingly, their physical and molecular properties, including their Soly IDs, chromosome information, genomic positions, lengths of both the CDSs and proteins, theoretical isoelectric points, and molecular weights were predicted (Table S2). All identified C2H2-ZF genes were found to encode proteins varying from 96 to 1,178 amino acids, including a few exceptionally longer (Solyc04g074250.2.1) or smaller proteins (Solyc00g015730.2.1). Correspondingly, their molecular weights varied from 11.021 kDa (Solyc00g015730.2.1) to 128.573 kDa (Solyc04g074250.2.1), and their predicted isoelectric points (pI) ranged from 4.49 (Solyc10g085560.2.1) to 10.62 (Solyc00g014800.1.1). Detailed information of tomato C2H2-ZF family is shown in Table S2.

Comparative phylogenetic analysis of the C2H2-ZFs

To explore the phylogenetic relationships of the C2H2-ZFPs between tomato and the model plant Arabidopsis, a NJ phylogenetic tree was firstly constructed from an alignment of the 104 tomato C2H2-ZFPs and 97 Arabidopsis C2H2-ZFPs (Table S3). The sequences were clustered into six major groups (Fig. 1A, Group 1-6, Table S3), each of which could be further subdivided into different subgroups according to their bootstrap values. Furthermore, the tomato C2H2-ZFPs were tightly grouped with the Arabidopsis C2H2-ZFPs (Groups 1 to 4, and Group 6; Fig. 1A), except for the Group 5 (Fig. 1A), indicating that we could assume the putative function of these genes in tomato based on the Arabidopsis C2H2-ZFPs in same group. Subsequently, we reconstructed the phylogenetic distribution of the 104 C2H2-ZF genes using different methods to clearly elucidate the relationships between the C2H2-ZF family genes in tomato. Results showed that the phylogenetic tree constructed by different methods displayed the differences, such as NJ tree (Fig. 1B), ML tree (Fig. S1A), UPGMA tree (Fig. S1B) and MP tree (Fig. S1C). For the NJ phylogenetic tree, 104 C2H2-ZFPs were classified into nine groups (I to IX; Fig. 1B). Among them, four genes in Group VII (Fig. 1B) were all grouped in Group 1 (Fig. 1A), some genes in Group VII, VIII and IX (Fig. 1B) was closely related to Group 1 and 2 (Fig. 1A), and Group II, III, IV (Fig. 1B) was closely related to Group 3 (Fig. 1A), while a small number of C2H2-ZFPs were clustered into the different groups in tomato (Figs. 1A and 1B). These results suggested that the same clades of phylogenetic tree were likely to represent the closest homologous gene pairs between tomato and Arabidopsis, and C2H2-ZFs have also undergone sequence diversification independently in different organisms. In addition, 104 C2H2-ZFs were irregularly distributed between these groups in tomato (Fig. 1B), such as Groups II, VIII, and IX were the main clades, containing 18, 26, and 28 genes, respectively, while Groups I, IV, and VII were the smallest clades, containing three, three, and four tomato C2H2-ZFs, respectively (Fig. 1B).

Conserved domain analysis of the C2H2-ZFPs

The ZF domain structure is X2-Cys-X(2-4)-Cys-X12-His-X(3-5)-His, where X represents any amino acid and the numbers indicate the consensus spacing between the conserved amino acid residues. This is highly conserved and essential for the ZF configuration and loop stability. To investigate the characteristics of the C2H2-ZF domains in tomato, the 104 C2H2-ZFPs were analyzed using the NCBI Conserved Domain Database (CDD; https://www.ncbi.nlm.nih.gov/cdd/). Multiple protein sequence alignments revealed that the ZF domains in all of these proteins were highly conserved. Among them, 83 C2H2-ZFPs (about 79.81%) contained a X2-Cys-X2-Cys-X12-His-X3-His sequence, seven C2H2-ZFPs (about 6.73%) had a X2-Cys-X2-Cys-X12-His-X4-His sequence, and the remains showed different types of C2H2-ZF domains (Fig. 2). In addition, 50 C2H2-ZFPs with a plant-specific conserved domain ‘QALGGH’ were identified in tomato and mainly classified into Groups VIII and IX (Fig. 2), which were fewer than that in Arobidopsis (64) (Englbrecht, Schoof & Bohm, 2004), rice (65) (Agarwal et al., 2007), poplar (62) (Liu et al., 2015) and foxtail millet (97) (Muthamilarasan et al., 2014), respectively.

Analysis of the exon–intron structures and conserved motifs in the tomato C2H2-ZFs

The Gene structural diversity within a gene family can be used as an evolutionary marker (Zhu et al., 2018). To gain further insights into the structural diversity of tomato C2H2-ZFs, we compared the exon–intron structure of each of the C2H2-ZFs. The number of introns in the C2H2-ZFs varied from 0 to 10. Based on the number of introns, 65 C2H2-ZFs (62.5%) was intronless, 29 C2H2-ZFs had 1 to 3 introns (27.9%), and the remaining C2H2-ZFs contained more than four introns (9.6%). In general, genes in the same group shared similar intron/exon arrangements in terms of intron numbers and exon length. For example, most C2H2-ZFs in Groups III, IV, VI, VII, VIII and IX genes had zero to one intron, and in Group II had two to three introns, which had the normal intron length (Fig. 3, Table S2). In contrast, the gene structure appeared to be more variable in groups I and V, which had striking distinctions in the exon/intron structure variants (Fig. 3, Table S2). In addition, we compared the exon–intron structure of each of the C2H2-ZFs in the same cluster of the different phylogenetic tree constructed by NJ (Fig. 3), ML (Fig. S1A), UPGMA (Fig. S1B) and MP (Fig. S1C), respectively. We found that genes in the same group had the highly similar intron/exon arrangements in the NJ phylogenetic tree (Fig. 3), indicating that the gene structure patterns were consistent with the NJ phylogenetic analysis.

Using the MEME tool, a total of 10 conserved motifs were identified in the C2H2-ZFPs, and the lengths of these conserved motifs varied from 7 to 31 amino acids (Fig. 4, Fig. S2). Among them, Motif 1 was widely detected in all C2H2-ZFPs, corresponding to the C-X2-C-X12-H-X3-H single ZF structure, which located at the N-terminal region of C2H2-ZFPs in groups II, and VI-IX (Fig. 4). Some groups also contained specific motifs; for example, Motifs 2 and 3 were present in the N-terminal region of Group II and the C-terminal region of many members in Groups III and IV; while Motif 4 was only detected in the N-terminal region of Group II, indicating that they may be relevant to the specific functions of these genes. The members of Groups I and V–IX have relatively simple motif patterns in comparison with Groups II to IV, implying the possible functional divergence of the C2H2-ZF genes in tomato.

Taken together, the structure and motif conservation within the C2H2-ZF genes supports the results of the NJ phylogenetic analysis. Variations in the motif compositions between subfamilies might be explained by their functional diversification.

Chromosomal locations and gene duplication events in the C2H2-ZFs

To gain insight into the organization of these genes in the tomato genome, the 104 C2H2-ZFs were mapped onto their respective chromosomes, which acquired from tomato genome database (The Tomato Genome, 2012). The C2H2-ZF genes were unevenly distributed throughout the tomato genome and generally more abundant at the both ends of the chromosomes (Fig. 5, Table S2). Chromosome 6 contained the largest number of C2H2-ZFs (16), followed by chromosomes 1, 2, 3, 4, 5, 9, 10, and 11, each of which contained 7 to 12 C2H2-ZF s. Chromosomes 0 (random chromosome) and 12 contained two and one genes, respectively. In addition, the duplication events of C2H2-ZFs were also analyzed in tomato genome since gene replication play an important role in genomic expansions and realignments. We identified 15 pairs of tandem-duplicated gene pairs (with two or more homologous genes within 100 kb region) located on chromosomes 1, 2, 4, 5, 9, 10, and 11 (Fig. 5), respectively. The closely related clustered sequences of Groups VI, VIII, and IX are mainly located on chromosomes 1, 2, 5, 9, 10, and 11, suggesting that the expansion of this gene family may have occurred via localized or intra-chromosomal duplication. In addition to tandem duplications, seven segmental duplication events were detected to scatter in eight chromosomes (Fig. 5). Furthermore, many homologous genes were located in different chromosomes in tomato, supporting the high conservation of the C2H2-ZF gene family.

Functional annotation of C2H2-ZFs

All 104 C2H2-ZFs were subjected to GO enrichment for investigating their functional annotation. As a result, a total of 43 C2H2-ZFs were mapped on the GO database, resulting in 285 annotations (Table S5), which were distributed across three ontology categories (biological processes, cellular components, and molecular functions) with 30 function terms (Fig. S5). Furthermore, most tomato C2H2-ZFs were annotated as being associated with meal ion binding, DNA-binding transcription factor activity, nucleus and regulation of transcription, DNA-templated, respectively (Fig. S5).

Expression patterns of the C2H2-ZFs in various tissues and organs

Using the TomExpress database (available online: http://tomexpress.toulouse.inra.fr/), we investigated the expression levels of the 104 C2H2-ZFs in various tissues of the tomato cultivar Micro Tom, including its roots, leaves, flower buds, petals, flowers, fruits, flesh, and peel, which were visualized using a heatmap (Fig. 6). The C2H2-ZFs had different expression patterns across the various tissues. Firstly, the clear differences in expression between the genes of Group II, III and V in comparison with those of the other groups were consistent with their different gene structures and the phylogenetic tree analysis (Figs. 3 and 6). In addition, we found that a small few genes typically had relatively high expression levels across various tissues and organs at different development stages, such as Solyc01g099340.3.1 in Group II, Solyc09g009030.3.1 in Group V, Solyc04g077980.1.1 in Group VIII, and Solyc06g075780.1.1 in Group IX, while most genes of Groups IV, VI, VII, and VIII had relatively low expression levels (Fig. 6). Some genes, including Solyc05g054030.3.1, Solyc05g055500.1.1, Solyc07g063970.3.1, Solyc09g011110.1.1, Solyc10g080600.2.1, and Solyc11g066400.1.1, were specially expressed in the different tissues, and genes (Solyc01g005130.2.1, Solyc04g014540.2.1, Solyc06g075780.1.1, and Solyc07g006880.1.1) in Group IX showed similar expression patterns during the fruit development (Fig. 6). However, the expression patterns of a few gene pairs, including Solyc03g025440.3.1, Solyc04g056320.2.1, Solyc06g065440.1.1, and Solyc11g017140.2.1, were significantly different in Groups I and III, although they were paralogous genes (Figs. 1B and 6).

Expression of the C2H2-ZFs under heat stress

To improve the reproducibility and reliability of the qRT-PCR analysis obtained in this study, three previously reported tomato internal reference genes, including SlEF1-α (Solyc06g005060), SlACT (Solyc03g078400), and SlUBI3 (Solyc01g056940), were primarily detected for the their quantification stability in the leaves, stems, and roots under the control and heat stress by qRT-PCR method. The results showed that SlEF1-α (Fig. S3A) displayed relatively greater stability than SlACT and SlUBI3 among the different tomato tissues (Figs. S3B and S3C, Table S4), which were used for further analysis in this study. Many of the C2H2-ZFs previously described in other crops were involved in responding to abiotic stresses, such as drought, heat, and salt stress (Muthamilarasan et al., 2014; Sakamoto et al., 2004; Sun et al., 2010; Wang et al., 2018). To confirm the potential roles of candidate tomato C2H2-ZFs in heat stress, the expression profiles of 34 C2H2-ZFs were finally detected in leaves, stems and roots of tomato under heat stress (42 °C), which might not be markedly induced by the day/night cycle, except the Solyc10g085560.2.1 (Fig. S4). Combined with the expression phylogenetic analysis, we found that genes showed the diversity expression patterns under heat stress (Figs. 7 and 8). Of them, the expression levels of 18 genes analyzed from Group II to VII were found to differ in roots, stems and leaves, under heat stress (Figs. 7A–7S). Of these, Solyc02g062940.3.1 and Solyc02g068990.3.1 were sustained up-regulated in the roots, stems, and leaves during heat stress treatment (Figs. 7F and 7R), and Solyc02g068980.1.1 was significantly up-regulated in the roots and stems with the similar expression patterns (Fig. 7P), but Solyc08g006470.3.1 was only significantly up-regulated by the heat treatment in the roots (Fig. 7L). For the largest classification groups VIII and IX, the expression levels of seven and nine genes from Group VIII and IX, were also investigated (Fig. 8). In Group VIII, Solyc05g012490.1.1 and Solyc04g005440.1.1 might be hard to induce by heat treatment with lower expression levels in leaves (Figs. 8A and 8B), while Solyc10g084910.2.1 and Solyc09g011110.1.1 showed the similar expression patterns and only significantly up-regulated in the roots (Figs. 8C and 8D). In contrast, Solyc06g053720.2.1 and Solyc12g088390.1.1 were significantly suppressed in the roots, and leaves under heat treatment (Figs. 8E and 8F). Notably, the heat treatment induced significantly high levels of Solyc01g107170.2.1 expression in all tissues (Fig. 8G). In Group IX, only Solyc01g005190.1.1 was significantly up-regulated in all tissues after heat treatment (Fig. 8M), while Solyc06g075780.1.1 and Solyc05g006310.2.1 were strongly up-regulated in the roots and stems (Figs. 8I and 8K). Although these genes showed various expression patterns under heat treatment, most of them in the identical group had also similar expression patterns in the same tissues, such as Group II genes (Solyc06g075250.3.1, Solyc03g121660.3.1, and Solyc02g062940.3.1) in stems (Figs. 7D, 7E and 7F) and (Solyc10g084180.2.1, Solyc06g075250.3.1, and Solyc02g062940.3.1) in leaves (Figs. 7C, 7D and 7F), Group VII genes (Solyc02g068980.1.1, Solyc02g068990.1.1, and Solyc08g059770.1.1) in roots and leaves (Figs. 7P, 7R and 7S), and Group VIII genes (Solyc05g012490.1.1 and Solyc04g005440.1.1) in the stems (Figs. 8A and 7B), which were significantly up-regulated under heat treatment, and Group VIII genes (Solyc05g012490.1.1 and Solyc04g005440.1.1, Solyc10g084910.2.1 and Solyc09g011110.1.1) in stems or leaves (Figs. 8A–8D) with lowly expression levels, while Group VIII genes (Solyc06g053720.2.1 and Solyc12g088390.1.1) were significantly suppressed in roots and leaves (Figs. 8E and 8F), and so on.

Taken together, qRT-PCR was performed for 34 selected genes from different cluster groups under heat stress, and showed the various expression patterns in the roots, stems, or leaves. Among them, seven genes were strongly up-regulated in the roots, stems, and leaves, while one gene were down-regulated; and five genes up-regulated both in roots and stems, and five genes specially expressed in the roots, and so on (Figs. 7 and 8). Hence, these results suggest that these C2H2-ZF genes might be associated with the heat stress during the seedlings development in tomato.