Identification and Validation of Selected Universal Stress Protein Domain Containing Drought-Responsive Genes in Pigeonpea (Cajanus cajan L.)

Pigeonpea is a resilient crop, which is relatively more drought tolerant than many other legume crops. To understand the molecular mechanisms of this unique feature of pigeonpea, 51 genes were selected using the Hidden Markov Models (HMM) those codes for proteins having close similarity to universal stress protein domain. Validation of these genes was conducted on three pigeonpea genotypes (ICPL 151, ICPL 8755, and ICPL 227) having different levels of drought tolerance. Gene expression analysis using qRT-PCR revealed 6, 8, and 18 genes to be ≥2-fold differentially expressed in ICPL 151, ICPL 8755, and ICPL 227, respectively. A total of 10 differentially expressed genes showed ≥2-fold up-regulation in the more drought tolerant genotype, which encoded four different classes of proteins. These include plant U-box protein (four genes), universal stress protein A-like protein (four genes), cation/H(+) antiporter protein (one gene) and an uncharacterized protein (one gene). Genes C.cajan_29830 and C.cajan_33874 belonging to uspA, were found significantly expressed in all the three genotypes with ≥2-fold expression variations. Expression profiling of these two genes on the four other legume crops revealed their specific role in pigeonpea. Therefore, these genes seem to be promising candidates for conferring drought tolerance specifically to pigeonpea.


INTRODUCTION
Abrupt climate changes and unavailability of sufficient water supply can severely affect the productivity of agriculturally important crops. Additionally, frequent exposure of environmental stresses such as drought is adversely affecting the plant growth and yield. Drought can occur at any stage of plant growth and the degree of yield loss depends on the onset time, intensity and duration of stress (Hu and Xiong, 2014). Pigeonpea is usually grown under marginal environments that are often subjected to water stress at different stages of growth and development. Even for short-duration varieties, yield gets affected due to water stress during late flowering and early pod development stages (Lopez et al., 1996). During seed hardening, the crop requires considerable amount of water and at this crucial stage unavailability of water often causes terminal drought. Despite having a deeper root system, drought is still one of the major yield-limiting factors, especially at critical seedling and reproductive stages of pigeonpea (Saxena, 2008). There has been a rousing progress made in developing drought-tolerant pigeonpea genotypes, but still it is difficult to meet the conditions arisen due to climate change. It is feasible to develop drought tolerant varieties through genomics-assisted breeding that would facilitate yield stability under water-deficient conditions (Varshney et al., 2014).
Since drought is a complex trait and is controlled by multigenes, identification of candidate genes and understanding the molecular mechanism associated with drought tolerance in pigeonpea is critical. Many studies have been carried out in model plants to identify candidate genes associated with drought response (see Mir et al., 2012). In pigeonpea, ample amount of genomics resources has been developed which can be deployed to identify candidate drought tolerant genes specific to pigeonpea. Moreover, the pigeonpea genome sequence reported 111 homologous sequences corresponding to universal drought-responsive protein sequences from the Viridiplantae . Similarly, the development of comprehensive transcriptome assembly (Kudapa et al., 2012) and the identification of genes involved in abiotic stresses tolerance have been reported Sekhar et al., 2010;Saxena et al., 2011;Deeplanaik et al., 2013).
Functional characterization of genes involved in different stress-responsive pathways such as photosynthesis and carbohydrate metabolism (Basu et al., 1999), related to stress-responsive transcription factors (Nakashima et al., 2009), signal transduction and regulatory compounds (Ramanjulu and Bartels, 2002;Sreenivasulu et al., 2007) gives an insight into the mechanisms adopted by plants to cope with drought stress. In this context, using bioinformatics approach, a total of 32 drought-responsive ESTs were retrieved from seven plant genera namely, Glycine, Hordeum, Manihot, Medicago, Oryza, Pinus, and Triticum (Isokpehi et al., 2011). Similarly, in soybean, 32 drought responsive genes involved in 17 metabolic pathways were identified and were validated in pigeonpea to know their association with drought stress (Deeplanaik et al., 2013).
The present study involves in-silico identification of selected universal stress protein domain containing drought-responsive genes. The qRT-PCR validation of these genes was carried out on pigeonpea genotypes with different levels of drought tolerance. Drought stress was imposed on all the selected genotypes and compared with well-watered controls to validate the candidate genes involved with drought tolerance in pigeonpea. The genes were also validated using a tolerant and a susceptible genotype each from four legumes namely, chickpea, groundnut, common bean, and cowpea. The identified candidate genes in future, can be functionally validated using transgenic approaches.
Additionally to utilize the identified drought tolerant genes, markers can be developed using haplotype analysis approach, which will accelerate crop yield even under drought stress conditions.

Drought Stress Treatment and Tissue Harvesting
Seeds were thoroughly washed with DEPC treated water, sown in 3 inches plastic pots (one seed per pot) filled with autoclaved black soil, sand, and vermicompost (10:10:1 v/v) mixture. All the plants were grown under controlled conditions in three replications. For imposing drought stress, slow drought (dry down) stress was imposed on the plants when they reached 22 days old seedling stage. A calculated amount of water was added to each pot, which was weighed at regular intervals. Control plants were maintained throughout at 80% relative water content (RWC) whereas stressed plants were dried down gradually to 20% RWC. The intensity of the drought stress was measured by recording the transpiration ratio (TR) on a daily basis. Stressed

RNA Isolation and cDNA Synthesis
Total RNA was isolated from all the frozen root samples using TRIzol (Invitrogen, USA) and was purified using DNase (Qiagen, GmbH, Germany) through an RNeasy Plant Mini kit according to the manufacturer's instructions. The concentration of total RNA was checked using NanoDropND-1000 (NanoDrop Technologies, USA) and RNA integrity was assessed on 1% denaturing formaldehyde agarose gel. cDNA was prepared using one microgram of total RNA using the SuperScript R III RT enzyme (Invitrogen, USA).

In-silico Identification of Drought Responsive Genes
To predict drought responsive genes in pigeonpea, gene set annotated with pigeonpea genome assembly v5  was downloaded from International Initiative for Pigeonpea Genomics (IIPG: http://www.icrisat.org/gt-bt/ iipg/Home.html). In addition, HMM profile of USP domain (PF00582) was retrieved from Pfam database (http://pfam. sanger.ac.uk/; Finn et al., 2010). The whole gene set was searched using "hmmsearch" program of HMMER 3.0 with HMM profile of the USP domain. Genes that detected sections encoding the USP domain, above the default inclusion threshold with statistically significant domain architecture were selected as USP domain-encoding genes. Sequences were also subjected to protein homology search using BLASTP against Swiss-Prot and TrEMBL databases to further determine the identity of the selected genes containing USP domains. The identities obtained from the two databases were used in UniProtKB database to retrieve the protein names, location, biological pathways, and gene ontology identity with the help of an in-house Perl script. Gene ontology enrichment analysis was performed using the BiNGO tool with p-value cutoff of ≤ 0.05 (Maere et al., 2005).

Primer Designing and qRT-PCR
Gene specific primer pairs were designed from the exonic regions of the selected genes. Primer3 software (http://probes. pw.usda.gov/cgi-bin/batchprimer3/batchprimer3.cgi) was used for primer designing using the following criteria: annealing temperature (Tm) in the range of 55-60 • C with an average of 57 • C, amplicon size of 150-200 bp, primer length of 20 ± 5 bp and GC% of 50 ± 5 (Supplementary Table 1). All the designed primer pairs were custom synthesized by MWG (MWG-Biotech AG, Bangalore, India).
qRT-PCR was carried out using ABI SYBR R GREEN PCR reaction on an ABI Fast7500 System (Applied Biosystems, Foster City, CA, USA) following manufacturer's instructions. PCR conditions maintained for all qRT-PCR reactions were 2 min at 50 • C, 10 min at 95 • C, and 40 cycles of 15 s at 95 • C and 1 min at 60 • C. The melt curve analysis was conducted for all 51 primer pairs. Only after confirming the observed single peak with all the selected tissue samples, primers were used further for qRT-PCR analysis.
Each reaction was carried out in three biological and two technical replicates along with no template control. The differential expression values of drought responsive genes were normalized with ACT1 as the reference gene (Supplementary Table 2). Statistical comparison between data obtained from three genotypes (ICPL 227, ICPL 8755 and ICPL 151) was performed using Tukey's post-hoc multiple comparison test using SPSS (version 16.0) whereas Student t-test was used to compare tolerant and susceptible genotypes of other four legume crops.

Identification of Drought Responsive Genes
The gene set consisting of 48,680 gene models  was used to search the USP domain-encoding genes. As a result, 71 genes were found to encode USP domain, of which 51 were identified to be above the inclusion threshold with E < 0.01 (Supplementary Table 3). Of these genes, 49 also showed identity to the 111 drought responsive genes reported earlier in pigeonpea .

Functional Classification of Drought Responsive Genes
To classify the 51 drought responsive genes based on their functional annotations, BLASTP search was performed against Swiss-Prot and TrEMBL databases. This analysis revealed that about 25.5% of the genes were classified as related to molecular functions such as catalytic activity (8%), transporter activity (5%), and binding (6%), whereas 27.4% of the genes were found to be related to cellular component such as ubiquitin ligase complex (5%), membrane (6%), organelle (5%), membrane part (5%), plastid part (1%), and cell part (11%). The genes involved in biological processes formed 47% and included response to stress (13%), metabolic process (8%), cellular process (12%), homeostatic process (3%), single-organism process (6%), localization (5%) and establishment of localization (5%). Gene ontology (GO) term enrichment performed using BiNGO for these 51 genes as visualized in Cytoscape has been presented in Figure 1. Detailed information of the corresponding protein name, GO term and ontology identities of these genes has also been provided in the Supplementary Table 4 and Supplementary Figure 1.
FIGURE 1 | Gene ontology (GO) enrichment of drought responsive genes. Genes showing differential regulation were analyzed using BiNGO and the biological process terms showing significant enrichment are presented. The colors shades represent the following significance level; white-no significance difference; yellow P = 0.05, orange P < 0.0000005.

Differentially Expressed Drought Responsive Genes
Gene specific primer pairs were designed from the exonic regions of 51 selected drought responsive genes for validation using qRT-PCR (Supplementary Table 3). Hierarchical cluster analysis of expression data of these genes showed a range of differential gene expression among three genotypes (ICPL 151, ICPL 8755, and ICPL 227) under stressed and controlled conditions. This analysis revealed that the less drought tolerant (LDT) genotypes, ICPL 151 and ICPL 8755 clustered separately from the more drought tolerant (MDT) genotype, ICPL 227 (Figure 2). Expression data was analyzed further in two different ways: (1) comparison between stressed and control samples for each genotype and (2) pair-wise comparison between genotypes. Expression analysis for each genotype with respect to stressed and control samples identified 18 genes in ICPL 227, six genes in ICPL 151 and eight genes in ICPL 8755 with significant expression variation. Genes with more than two-fold expression difference in each of these genotypes have been listed in Table 2. However, the expression analysis among the genotype pairs, ICPL 227 with ICPL 8755 (Supplementary Table 5) and ICPL 227 with ICPL 151 (Supplementary Table 6) has identified 11 genes in each case.

Comparative Expression Profiling of Candidate Genes across Legumes
Among the 10 common differentially expressed genes among MDT and LDT genotypes, two genes namely, C.cajan_29830 and C.cajan_33874 showed marked up-regulation in all the three drought tolerant genotypes (Supplementary Figure 2). These two genes were further selected for validation using tolerant and susceptible genotypes in chickpea, groundnut, common bean, and cowpea. For gene normalization, glyceraldehyde 3-phosphate dehydrogenase (GAPDH) for chickpea (Garg et al., 2010), alcohol dehydrogenase (ADH) for groundnut (Reddy et al., 2013) while β-tubulin for common bean and cowpea (Eticha et al., 2010) were used as internal control (Supplementary Table 2 and Supplementary Figure 3).
To perform normalization, the reference genes which were reported to be stable for each of the legume crops were considered. In terms of expression profiling, in chickpea, ICC 4958 (tolerant genotype), showed 0.53, 0.14, and ICC 1882 (susceptible genotype), showed 0.23-, 1.01-fold differential gene expressions for C.cajan_29830 (Figure 4) and C.cajan_33874 (Figure 5), respectively. Similarly in groundnut, these two genes showed an expression variation of 0 Overall, the gene expression variation between the tolerant and susceptible genotypes of the four legumes for the selected genes was less than 2-folds. Thus, the present study implies that these two genes might specifically be involved in contributing drought tolerance in pigeonpea.

DISCUSSION
This study has utilized genome sequence information for selecting genes encoding proteins containing USP domain and were validated for their role in drought tolerance in pigeonpea using qRT-PCR. Genes encoding protein with USP domain are known to be involved in a myriad of stress responses and any mutation in these genes may cause loss of efficacy in combating stresses (Drumm et al., 2009;Isokpehi et al., 2011;Shokry et al., 2014). USP domain has been found to be evolutionary conserved in a number of crop species such as cassava, soybean, finger millet, and peanut (Govind et al., 2009;Deeplanaik et al., 2013;Turyagyenda et al., 2013). Earlier studies have provided evidences that stress-responsive genes encoding proteins with USP domain are useful in stress signal perception and subsequently lead to functionally efficient proteins. These proteins have been found to be involved in protecting cellular structures and cell molecules under stress conditions (Waditee et al., 2002;Majee et al., 2004;Dastidar et al., 2006;Govind et al., 2009). For instance, under water deficit conditions, out of 50 genes selected in peanut, only HSP70 gene showed association with drought stress response (Govind et al., 2009). Similarly, 10 genes conferring drought tolerance were characterized in cassava (Turyagyenda et al., 2013). In the FIGURE 4 | Expression analysis of identified drought responsive gene, C.cajan_29830 in the root tissues of tolerant and susceptible genotypes of the four legumes. The expression analysis was carried out using chickpea genotypes, ICC 4958 (tolerant) and ICC 1882 (susceptible), CSMG 84-1 (tolerant) and ICGS-76 (susceptible) for groundnut, IT93K503-1 (tolerant) and UC-C B46 (susceptible) for cowpea, and BAT 477 (tolerant) and DOR 364 (susceptible) for common bean, in comparison to the pigeonpea genotypes, ICPL 151 (LDT genotype) and ICPL 227 (MDT genotype).
case of pigeonpea, homology search provided 71 genes encoding USP domain, of which 51 genes with pure domain architecture were selected for further validation in pigeonpea. During water stress conditions, plant responds at both cellular as well as molecular level by accumulating osmolytes and proteins involved in stress response and/or tolerance (Yamaguchi-Shinozaki and Shinozaki, 2006). Stress response at cellular level such as cell proliferation, differentiation, stomatal closure, repression of cell growth generally lead to induced expression of drought responsive genes (Yamaguchi-Shinozaki and Shinozaki, 2005). Among the selected 51 genes, majority of the genes were found to be related to response to stress followed by cellular processes.
The three pigeonpea genotypes selected for expression profiling in the present study exhibited varying degree of tolerance to drought stress (Lopez et al., 1996;Saxena et al., 2011). The expression variation of candidate genes in stressed tissues can be compared with well water controls at specific time frame (VanGuilder et al., 2008) using qRT-PCR. A set of 10 genes, which were identified in the MDT and LDT genotypes showed large difference in the relative transcript abundance. These genes were found to be related to plant U-box proteins, cation/H(+) antiporter (CHX) proteins, uncharacterized protein and universal stress protein A-(uspA) like protein.
Plant U-box E3 ligases have been found to be involved in enhanced drought, salinity, cold and heat tolerance in Arabidopsis thaliana (Lyzenga and Stone, 2011). Whereas, Ubiquitin-protein ligases (E3s) determine the substrate specificity of ubiquitylation and plays an important role in protein post-translational modification in higher plants (Liu and Walters, 2010). Based on the structure, Ubiquitin-protein ligases (E3s) has been classified into two families, the HECT and RING-finger (U-Box) families (Hatakeyama and Nakayama, 2003). Molecular and cellular characterization of U-Box proteincoding genes in hot pepper (Cho et al., 2006) and Arabidopsis (Cho et al., 2008) has clearly demonstrated the role of U-Box protein-coding genes in drought tolerance. The four plant U-Box protein-coding genes (C.cajan_26230, C.cajan_39705, C.cajan_09181, and C.cajan_30211) identified between the MDT and LDT genotypes showed significant differences in the relative transcript abundance.
One gene (C.cajan_46779) belonging to cation/H(+) exchanger (CHX) group was also found to be differentially expressed between the two genotypes studied. In Arabidopsis, CHX gene family was found to play an important role in osmotic adjustment and K+ homeostasis (Sze et al., 2004). Therefore, the finding suggests that the gene, C.cajan_46779 belonging to CHX gene family may also be playing an important role during drought stress condition in pigeonpea. Another gene, C.cajan_08737, annotated as uncharacterized protein, was also found to be differentially expressed among MDT and LDT genotypes suggesting its role in drought tolerance. Another class of genes (C.cajan_13768, C.cajan_23080, C.cajan_29830, and C.cajan_33874) having differences in the relative transcript abundance between MDT and LDT genotypes was found to have the uspA domain. The uspA domain is also known to play a vital role in survival during cellular growth arrest. Genes belonging to this domain help in oxidative stress resistance and initiates defense against superoxide-generating agents (Nachin et al., 2005).
Pigeonpea is one of the most drought tolerant legume crops (Varshney et al., 2009). Therefore, the expression variation of the two differentially expressed genes, C.cajan_29830 and C.cajan_33874 identified across all the three pigeonpea genotypes was studied across four other legumes. For these two genes (C.cajan_29830 and C.cajan_33874), chickpea (78.5 and 82%), groundnut (77.6 and 79.5%) and common bean (85.9 and 86.4%) genes showed high sequence identity with pigeonpea, respectively. Based on this observation, two genes encoding Universal stress protein A-like protein seem to be highly conserved among the legumes studied. However, in the four legume crops, namely chickpea, groundnut, common bean, and cowpea, these two genes did not show any expression variation between the drought tolerant and susceptible genotypes. This may also be due to different strategies acquired by different legumes for drought adaptation mechanism (Mir et al., 2012). This study also showed the involvement of U-Box protein-coding in having some specific role in drought tolerance mechanisms in pigeonpea. Many such conserved U-Box protein-coding genes were found to be over-expressed in different plant species to enhance drought tolerance (Lyzenga and Stone, 2011).
Thus, expression analysis of the 51 drought responsive genes has provided a set of 10 genes belongs to plant U-Box proteins, cation/H(+) antiporter (CHX) proteins, uncharacterized proteins and universal stress protein A-(uspA) like protein. This candidate gene-based approach could provide useful insights into the molecular mechanisms involved in drought tolerance in pigeonpea. Moreover, the identified genes can also be validated at sequence level in different genetic backgrounds to detect the presence of sequence variations for the development of gene-based marker(s) for crop improvement and development of more tolerant breeding lines/hybrids through genomics-assisted breeding.

ACKNOWLEDGMENTS
The authors thank United States Agency for International Development (USAID) for financial support for the research work. This work has been undertaken as part of the CGIAR Research Program on Grain Legumes. ICRISAT is a member of CGIAR Consortium.
Supplementary Figure 2 | Expression variation of two candidate genes (C.cajan_29830 and C.cajan_33874) between MDT and LDT genotypes. Differentially expressed genes were identified with ≥ 2-fold expression variation across the three pigeonpea genotypes, namely ICPL 227 (MDT genotype), ICPL 151, and ICPL 8755 (LDT genotypes). The different letters above the bars were considered as statistically significant between each other.