Computational deciphering of biotic stress associated genes in tomato (Solanum lycopersicum).

Tomato (Solanum lycopersicum) is one of the major vegetable plant and a model system for fruit development. Its global importance is due to its lycopene pigment which has anti-oxidative and anti-cancerous properties. Though > 1.5 M biotic stress associated ESTs of tomato are available but cumulative analysis to predict genes is warranted. Availability of whole genome de novo assembly can advantageously be used to map them over different chromosome. Further, available 0.14 M catalogued markers can be used to introgress specific desirable genes in varietal improvement program. We report here 57 novel genes associated with biotic stress of tomato along with 50 genes having physical location over different chromosomes. We also report 52 cis-regulating elements and 69 putative miRNAs which are involved in regulation of these biotic stresses associated genes. These putative candidate genes associated with biotic stress can be used in molecular breeding in the endeavor of tomato productivity along with its sustainable germplasm management.

Introduction

Plants have developed gradually to live in an environment where they are generally exposed to different types of stresses in combination. Being sessile, they have developed various mechanisms which allow them to detect very minute changes and respond to those stress conditions thus minimizing damage [1]. These mechanisms activate expression of various genes expressions thus enabling them to maximize the chances of their survival [2]. Prolonged exposure of plants to biotic and abiotic stresses lead to reduction in fitness and ultimately in productivity [3]. Biotic stress induces a strong pressure on plants due to the attack of pathogens [4]. The attacking pathogens can be bacteria, fungus, Viruses or Nematodes. Bacteria enter the cell via gas or water pores (stomata and hydathodes) or any kind of wounds whereas; Fungi can enter the plant cell via extended hyphae on top or in between the plant cells. Nematodes and aphids obtain its feed by inserting a stylet directly into a plant cell [5]. Plants recognize the invaders by sensing evolutionary conserved microbial molecular signatures, known as pathogen-associated molecular patterns (PAMPs) or microbe-associated molecular patterns (MAMPs) by plant pattern recognition receptors (PRRs) [6], [7], [8]. After identifying PAMPs, PRRs triggers an immune response known as PAMP-triggered immunity (PTI), which imparts protection against non-host pathogens and limits disease caused by virulent pathogen [9]. Pathogen of the host plants release effector proteins which targets PTI components for inhibiting plant defense [10], [11], [12], [13], [14]. These effector proteins are recognized by Resistance proteins released by resistance genes (R-gene). Thus plants respond to these effector proteins by effector triggered immunity (ETI) which is highly specific and is assisted by hypersensitive response (HR) and systemic acquired resistance (SAR) [15], [16].

Many major plants, including tomato (viz. Solanum lycopersicum), have been targeted by various pathogens which responds to their attacks by ETI [17]. Tomato, which belongs to Solanaceae family, is among the major vegetable plants and is the model system for fruit development [18]. It is a short lived herb and a protective food as it contains important nutrients such as lycopene, beta-carotene, flavonoids, Vitamins A, B & C and hydroxycinnamic acid derivatives [19], [20], [21]. From last few years, it has attained significant popularity because of its lycopene's anti-oxidative and anti-cancerous properties [22].

Tomatoes can be infected by around 200 diseases, which causes reduction in its production [23]. Many resistant genes (R-genes) in this plant has been revealed, whose proteins identify the avirulent proteins of pathogens and initiate the defense mechanisms. Some of the identified R-genes are Cf for Cladosporium fulvum, Ve for Verticillium dahlia and Cmm for Clavibacter michiganensis [24]. For preventing diseases in susceptible varieties various chemical, physical and biological methods are followed by the farmers to prevent the disease [25]. Since these methods cost high, are time consuming and causes health hazards so there raised the need for genetically modified crops. Thus, for the species prone to diseases, transgenic varieties have been developed. Be2 R-gene from pepper (which is a close relative of tomato) has been introduced into tomato which has been effective in providing resistance against bacterial spot disease [26].

Though > 1.5 million ESTs of tomato are available through various publications, but cumulative analysis to predict evidently existing gene model along with its number are yet to be attempted. Since whole genome based de novo assembly of tomato is available, thus these EST based gene prediction can be more realistic than genome based. Biotic stress resistant trait is one of the most important trait in tomato for productivity, thus its associated genes must be mapped over different chromosomes of tomato. Physical location of such genes can be of great advantage as > 0.14 M markers [27] are already catalogued chromosome-wise. Selection of flanking region markers can facilitate the introgression of specific desirable genes in varietal improvement program.

Materials and methods

ESTs collection and its assembly

In our study, firstly ESTs of various diseases of tomatoes were collected from NCBI [28]. These ESTs were for several bacterial, fungal, nematode and viral diseases. The gathered sequences were further checked for duplication by CD-HIT Software. CD-HIT compares all the sequences, separate the duplicates and cluster them into groups [29]. In CD-HIT 0.95 (95% identity) was used as clustering threshold and word size was taken as 5. The rectified sequences were then assembled into contigs using an assembly program EGassembler and Velvet programs. EGassembler is a programmed and user-customized analysis tool for cleaning, repeat masking, vector trimming, organelle masking, clustering and assembling of ESTs and other genomic fragments [30]. Velvet is a package that deals with de novo genome assembly and short read sequencing alignments [31]. For EGassembler minimum percent identity for an alignment was taken as 96% and overlap percent identity cutoff as 80. As far as Velvet is concerned, minimum contig length was settled at 100 and K-mer length was taken as 21.

Functional annotation and gene ontology analysis of contigs

The assembled contigs were further functionally annotated. Thus functional characterization and gene ontology study was carried out by using Blast2GO 2.8 [32]. In Blast2GO, Blastx program, with e-value cut off as 1.0E − 3 and number of blast hits as 20, was used to search against the protein database (Refseq Protein Database) for all six translated reading frames of generated contigs for each disease separately [33]. Mapping and Interpro functions of Blast2GO were used to describe the exact GO terms associated with annotations obtained from BLAST result. Result from blastx exhibited, the homologous sequences present in protein database corresponding to the contigs. These sequences were further categorized into gene ontology categories viz. Molecular Function, Biological Process and cellular components. Contigs with zero hits and no functional information were further carried for gene prediction each.

Gene prediction

Unannotated contigs from Blast2GO were further loaded in FGENESH program of Molquest software version 2.4.5 [34]. Molquest is a software package of various gene prediction programs. Genes were predicted using reference genome of Solanum lycopersicum and Nicotiana tabacum. As, Solanum lycopersicum and Nicotiana tabacum belong to the same family Solanaceae [35], and in Molquest (with version 2.4.5), genomic information of only these two plants of that family is present, so Nicotiana tabacum was also used as reference so that maximum genes of that family can be mined.

Promoter analysis of novel candidate genes

PLant Cis regulatory DNA Elements (PLACE) database [36] was used to validate and identify regulatory regions like cis regulatory elements, transcription factor, promoter regions in all the candidate genes predicted from FGENESH program.

Chromosome mapping

For the predicted genes, there locations on different chromosomes of Solanum lycopersicum, were identified with the help of Blast. Thereafter, these genes were mapped with the help of MapViewer available at NCBI (http://www.ncbi.nlm.nih.gov/mapview/).

MicroRNA prediction

MicroRNAs (miRNAs) are small RNA species that play important regulatory roles in various biological processes in plants. miRNAs are known to regulate the expressions of many stress-related genes [37]. miRNAs were obtained from miRBase (release 21) database and were used for identifying their target genes among predicted genes. miRBase database is an interface for annotated miRNA sequence data [38]. For target identification, miRANDA software was used with energy value as − 4 kJ/mol as cut off. The miRanda algorithm is based on a comparison of miRNAs complementarity to 3′UTR regions [39].

Results and discussion

The ESTs of several diseases affecting tomatoes, obtained from NCBI, were categorized broadly Bacterial, Fungal, Nematodes and Viral. Total 8581, 26,277, 41,832 and 1,441,496 were collected for bacterial, fungal, nematodes and viral diseases respectively. After removing duplicates number of ESTs reduced to 7101 for bacterial, 6611 for fungal, 21,541 for nematodes and 933,034 for viral diseases. Assembly was carried out with both EGassembler as well as Velvet assembler. Assembled contigs with EGassembler and Velvet for all the diseases are given in Table 1. It was observed that contigs assembled with EGassembler were larger in number as compared to those assembled with Velvet.

Table 1

ESTs and assembled contigs.

Disease	ESTs/nucleotide collected	ESTs after repeat masking	Contigs generated (EGassembler)	Contigs generated (Velvet)
Bacterial	8581	7101	333	23
Fungal	26,277	6611	119	8
Nematode	41,832	21,541	3874	183
Viral	1,441,496	933,034	116,149	9371

Contigs obtained from Velvet and EGassembler were further submitted for functional annotation. After completing blast, it was perceived that Contigs obtained from EGassembler had 120,136 hits, while from Velvet had 4382 hits. So, Contigs of Velvet assembler and EGassembler without blast hits were 5200 and 339 respectively. Disease wise details for with and without blast hits (by two different tools viz., EGassembler and Velvet) are given in Table 2. Thereafter mapping and annotation was done for the contigs for which blast hits were obtained. These annotated sequences were categorized into three types — cellular components, biological processes and metabolic functions. For Velvet on an average 62.25% of the contigs were participating in biological processes, 67% for metabolic functions and 35.25% as cellular components. Whereas in case of EGassembler Contigs 66.25% were involved in biological processes, 64.5% for metabolic functions and 47.5% for cellular components on an average. A detail of this functional characterization is given in Table 3. Various pie charts for functional characterization of the annotated contigs are given Fig. 1.

Table 2

Disease wise details for with and without blast hits.

Disease	EGassembler		Velvet
	Contigs without blast hits	Contigs with blast hits	Contigs without blast hits	Contigs with blast hits
Bacterial	10	323	3	20
Fungal	5	114	0	8
Nematode	209	3665	38	145
Viral	115	116,034	5159	4212

Table 3

Functional characterization of the annotated contigs.

Disease	EGassembler			Velvet
	Biological processes(in %)	Molecular functions(in %)	Cellular components(in %)	Biological processes(in %)	Molecular functions(in %)	Cellular components(in %)
Bacterial	71	68	66	90	95	45
Fungal	63	60	45	75	88	50
Nematode	62	61	38	60	61	34
Viral	69	69	41	24	24	12

Fig. 1

Functional characterization of the annotated contigs (represented as pie charts).

Gene prediction and novel gene identification

The unannotated Contigs were used as input in FGENESH program of Molquest. It was observed that when Solanum lycopersicum was used as reference genome total 58 genes were predicted for Velvet Contigs and 38 genes for Contigs of EGassembler. But when Nicotiana tabacum (most phylogenetically related species to Solanum lycopersicum present in Molquest) was used as reference genome then 88 genes were predicted for Contigs of Velvet and 48 genes for EGassembler Contigs. After comparing the predicted genes, 50 novel genes were identified for viral disease, were there and for nematode there were 14 genes and for fungal disease there was one novel gene and for bacterial diseases there was one novel genes (with two isoforms) were identified.

Promoter analysis was performed to identify the probable cis-acting DNA sequences which may be liable in regulating candidate gene expression. The cis-acting sequences are actually a part of the gene and act as the regulatory sequences which influence the expression of the gene which contains them. Although, these sequences mostly found just upstream of the TSS but they can also be present much further upstream, or on the 3′ end of the gene, or even within the introns and exons of a gene [40]. Thus, the cis regulatory elements in the promoter region of the candidate genes were acquired from PLACE which is a database of motifs found in plant cis-acting regulatory DNA elements from all previously published reports. For the novel genes, cis regulating elements that regulate the biological process in diseased conditions were identified. To analyze this, PLACE database was used with all novel genes sequences as query sequences. Analysis of cis-regulatory elements describes the candidate gene expression and corresponding functional transcription factor. It was found that many diseases share the common transcription factors during pathogen attack. Cis-regulatory motifs and functioning transcription factor for each predicted candidate gene is given Supplementary sheet I. Among all the candidate gene, transcription factors like DOF, bZIP, WRKY, RAV, MYC and MYB expressed majorly. Details of transcription factors for genes predicted are given in Supplementary sheet I.

Chromosome mapping helps in locating the genes of interest on the genome of an organism. When the chromosome mapping of predicted genes, done with the help of Blast and Mapviewer, the map shown in Fig. 2 was obtained (detailed locations of the all the genes are given in Supplementary sheet II. It was observed that out of 55 genes 12 were located on chromosome number 1. Chromosome numbers 4 and 11 possessed equal no of genes i.e. 8. Details of no. of genes on each chromosome are given in Table 4.

Fig. 2

Chromosome mapping of the predicted genes.

Table 4

No. of genes on each chromosome.

Chromosome no.	No of genes
ch01	12
ch02	0
ch03	3
ch04	8
ch05	5
ch06	5
ch07	1
ch08	0
ch09	5
ch10	0
ch11	8
ch12	3

For all the predicted genes, microRNAs were predicted using miRANDA and miRBase data. It was found that 6 miRNAs were predicted for bacterial genes and 3 miRNAs for fungal genes. For genes related to nematodes had 24 miRNAs while that of genes related to viral diseases had 60 miRNAs. Details of the all the miRNAs for every gene is given in Supplementary sheets III. For total 50 genes, 69 unique miRNAs were identified.

Discussion

Development of stress tolerant varieties is one of the biggest challenges in plant breeding. For this the pathways related to various biotic stresses need to be deciphered and for this interpretation the genes should be known, so that the connectivity of these genes can be cracked. So keeping in view this point, new genes were identified from the available EST data at NCBI. These ESTs were assembled into contigs and these contigs were then annotated. From the annotation, it was observed that maximum number of contigs played role in metabolic processes and cellular process for bacterial, fungal, nematode and viral contigs.

Further, from unannotated contigs, 50 novel genes were predicted. For these genes promoter analysis was carried out. Many Transcription factors were identified for all the genes. The highly expressed transcription factors correspond to their role in biotic stress. It has found in previous studies that Dof transcription factor, with highly conserved Dof domain, play major role in two biotic stresses (watermelon mosaic virus and downy mildew) in cucumber [41]. Studies describe that Dof is a plant specific transcription factor and contains conserved C2-C2 zinc finger which help it in binding DNA. Thus it is known as DNA-binding One Zinc Finger [42]. bZIP transcription factor, also known as Basic Leucine Zipper, is present in all eukaryotes and is known to regulate processes corresponding to pathogen attack in many crops like Arabidopsis, cotton and maize [43]. Among the highly expressed WRKY transcription factor is participating in many stress (whether biotic or abiotic) signaling pathways. Thus, it acts as a good candidate for various stress tolerance mechanisms [44]. MYB transcription factor represents a family of proteins that comprise a conserved domain, the MYB DNA-binding domain. It has been reported that in Arabidopsis, MYB encodes an activator of the hypersensitive cell death program in response to pathogen attack [45]. Thus presence of various cis regulatory elements on the promoter region of the novel genes indicates their possible involvement in biotic stress pathways in Tomato.

To find the location of the novel genes, chromosome mapping was done. It was observed that all the predicted genes related to viral pathogen are located on chromosome no 1 so it can be said that chromosome no 1 is majorly involved in pathways related to attack by viral pathogens in plants. Later, using the gene data, when miRNAs were predicted it was observed that 69 unique miRNAs were there which were involved in various processes in tomatoes.

Conclusion

Tomato is a very significant fruit crop and is being widely used due to its good nutritive value. It is affected by various pathogens as a result its production goes down. So being an important plant, there is a need to develop new resistant varieties which can survive under the stress condition. So using the genomics approach, with ESTs as an input data, various novel genes have been deciphered. The advance biotechnological methods need to be revolutionized to the new predicted genes for developing new resistant breeds. Cis regulatory elements and transcription factors study provides good insight of their role in corresponding stress condition, whose validation is further warranted.

The following are the supplementary data related to this article.

Supplementary sheet I

Cis regulatory elements of the predicted genes.

Supplementary sheet II

Details of locations of predicted genes over different Chromosomes.

Supplementary sheet III

Various miRNAs of the predicted genes.