Genome-wide profiling reveals extensive alterations in Pseudomonas putida-mediated miRNAs expression during drought stress in chickpea (Cicer arietinum L.)

Environmental and Experimental Botany 157 (2019) 217–227

Contents lists available at ScienceDirect

Environmental and Experimental Botany journal homepage: www.elsevier.com/locate/envexpbot

Genome-wide profiling reveals extensive alterations in Pseudomonas putidamediated miRNAs expression during drought stress in chickpea (Cicer arietinum L.) Ram Jatana,b, Shalini Tiwaria, Mehar H. Asifa,b, Charu Lataa,b, a b

T

⁎

CSIR-National Botanical Research Institute, Rana Pratap Marg, Lucknow, 226001, India Academy of Scientific and Innovative Research (AcSIR), Anusandhan Bhawan, 2 Rafi Marg, New Delhi, 110 001, India

A R T I C LE I N FO

A B S T R A C T

Keywords: Drought miRNA PGPR Regulation Target gene Stem-loop quantitative real- time PCR

Pseudomonas putida strain MTCC5279 (RA) is a plant growth promoting rhizobacteria which improves growth and development of plants by colonizing their root surface and also confers drought stress tolerance. Further, along with several stress responsive genes, microRNAs are also well known for their crucial involvement in drought stress response. Considering the above, in this study, we performed identification and expression profiling of small RNAs in chickpea roots either uninoculated or inoculated with P. putida RA under drought stress through high-throughput sequencing. We identified 923 conserved and 216 novel miRNAs from all four libraries. Out of these, 431 and 332 conserved miRNAs; and 60 and 39 novel miRNAs were found to be differentially expressed in RA and drought with RA-inoculated plants, respectively. A total of 4647 and 551 target genes involved in various biological, cellular and molecular processes were predicted to be targets of conserved and novel miRNAs, respectively. Stem-loop quantitative real-time PCR analysis of selected nine miRNAs showed their expression patterns in consistence with those obtained from Illumina-sequencing. Besides, the expression patterns of miRNAs showed inverse correlation to that of their respective targets. Our results indicate that RA plays a very crucial role in regulating the expression of miRNAs and their targets during drought stress in chickpea paving the way for further characterization and utilization of miRNAs in crop improvement programmes.

1. Introduction Chickpea (Cicer arietinum L.) is one of the most important cultivated legumes that serve as a key source of essential nutrients and dietary protein for millions of people globally (Varshney et al., 2013; Srivastava et al., 2015). It helps in maintaining the soil nitrogen levels and crop diversification due to its ability to fix atmospheric nitrogen through symbiosis with soil bacteria (Varshney et al., 2013; Srivastava et al., 2015). It is mostly cultivated in arid and semi-arid areas and few chickpea genotypes are also capable to acclimatise to such conditions, however, its production and productivity are reportedly affected severely by several environmental stresses (Talebi et al., 2013; Kohli et al., 2014). Drought stress significantly reduces chickpea yield by 7–51% (Khodadadi, 2013; Hajyzadeh et al., 2015). It is however, well known that plants have developed several morpho-physiological, biochemical and molecular mechanisms that lead to reprogramming of the expression of various genes/metabolites and small non-coding

⁎

regulatory RNAs to avoid or tolerate drought conditions (Ahuja et al., 2010; Lata et al., 2015; Ma et al., 2015; Shivhare and Lata, 2017; Tiwari et al., 2017). Various studies have reported that small RNAs play very critical roles in regulating the expression of numerous stress-responsive genes at both transcriptional and post-transcriptional levels in various organisms (Trindade et al., 2010; Sunkar et al., 2012). MicroRNAs (miRNAs) are 20 to 24 nt long small RNAs that regulate target gene expression in plants and animals (Voinnet, 2009; Trindade et al., 2010; Sunkar et al., 2012). miRNA-mediated post-transcriptional silencing of genes takes place by the cleavage of mRNA or translational repression, based on the degree of complementarity of the miRNAs with its cognate target mRNAs (Llave et al., 2002; Bartel, 2004; Voinnet, 2009; Sunkar et al., 2012). Functional analyses have demonstrated that various cellular, biological and molecular mechanisms were regulated by miRNAs in plants (Zhang, 2015; Xia et al., 2016). In addition, they have also been shown to play crucial roles in various abiotic stress responses

Corresponding author at: CSIR-National Botanical Research Institute, Rana Pratap Marg, Lucknow, 226001, India. E-mail address: [email protected] (C. Lata).

https://doi.org/10.1016/j.envexpbot.2018.10.003 Received 23 June 2018; Received in revised form 19 September 2018; Accepted 4 October 2018 Available online 05 October 2018 0098-8472/ © 2018 Elsevier B.V. All rights reserved.

Environmental and Experimental Botany 157 (2019) 217–227

R. Jatan et al.

Only unique reads after removing redundant reads were selected for further analysis. After that, reads from different libraries which mapped to plant t/rRNAs were discarded by using the “filter tool” of sRNA workbench and the remaining sRNAs reads of 18–32 nucleotides were used for identification of both conserved and novel miRNAs.

including drought (Ferdous et al., 2017; Wu et al., 2017; Aravind et al., 2017; Liu et al., 2017), salt (Macovei and Tuteja, 2012; Li et al., 2013; Yin et al., 2017), heat (Liu et al., 2015), cold (Barakat et al., 2012; Yang et al., 2017) and CO2 (Xue et al., 2018) in different crop plants. Recently, a total of 224 conserved and 60 novel miRNAs were identified from the root apex of chickpea under drought and salt stresses (Khandal et al., 2017). A previous genome-wide miRNA profiling in chickpea reported a total of 440 conserved and 178 novel miRNAs from seven diverse tissues/organs (Jain et al., 2014). Further, Pseudomonas putida strain MTCC5279 (RA), a root-colonizing PGPR has the ability to improve growth and development of plants (Srivastava et al., 2012) and ameliorate drought stress in chickpea (Tiwari et al., 2016; Jatan et al., 2018). However, until now, there has been no report towards drought stress-responsive miRNAs discovery in chickpea or in any other plant in the presence of a plant growth promoting rhizobacteria (PGPR) except one reported by our group where we have studied the expression of at least nine miRNAs and their respective target genes in chickpea upon RA inoculation under drought and salt stresses indicating the crucial role of RA in stress response (Jatan et al., 2018). Genome-wide identification and characterization of miRNAs were, however, still needed to have a wider perspective on their role in PGPR-facilitated drought stress response and tolerance in chickpea. Therefore, in this study, we have performed genome-wide small RNA sequencing for identification of conserved and novel miRNAs and determination of their expression patterns in the presence and absence of RA. Further, we have also analysed the expression profiles of the targets of selected miRNAs. The outcomes will facilitate the elucidation of the functions of miRNAs in RA-mediated drought stress adaptation mechanisms and establish a framework for the increased understanding of their role in drought stress alleviation in chickpea.

2.3. Identification and differential expression of conserved and novel miRNAs in chickpea To predict conserved and novel miRNAs, miRProf and miRCat (miRNA Categorization) pipeline of UEA sRNA workbench were used. The filtered reads which mapped with miRNA sequences available in miRBase (release 21; Griffiths-Jones et al., 2008) entries of other plant species using Bowtie alignment tool (Li and Durbin, 2009) were designated as conserved miRNAs. The miRCat was run with the default parameters to identify novel miRNAs. The location of novel miRNAs on the genome of chickpea was determined by using MapChart based on their genomic coordinates (Voorrips, 2002). RNA/folding annotation tool (Vienna RNA software package) of UEA sRNA workbench was used to establish the secondary structure of precursor sequences from the 100 nt flanking regions of aligned reads of the genome (Meyers et al., 2008). The minimum free energy (MFE) of precursor sequences was used to determine the stability of RNA secondary structure (Llave et al., 2002). To identify differentially expressed miRNAs, edgeR method was used (Robinson et al., 2010). The miRNA read counts of treated libraries and that of control library were used to calculate the fold change of miRNAs. miRNAs having log2 fold change ≤ ± 1.0 and pvalue < 0.05 were considered as significantly differentially expressed. 2.4. Prediction and functional annotation of the target genes The psRNATarget server was used to predict the potential target genes of the novel and conserved chickpea miRNAs with default parameters (Dai and Zhao, 2011). The putative functions of predicted target genes were identified by using chickpea genome annotation (Jain et al., 2013). The gene ontology (GO) enrichment analysis was performed by using PlantRefSeq and UniProt database for all predicted target genes.

2. Materials and methods 2.1. Plant materials and stress treatment A drought-tolerant chickpea cv. BG-362 was used for identification and expression analysis of miRNAs in the presence or absence of RA which was obtained from our lab repository. Chickpea seeds were surface sterilized and grown in hydroponics in a controlled growth chamber with and without 1% RA inoculation as described previously (Jatan et al., 2018). After 24 h of inoculation, trays designated for drought and drought with RA were supplemented with 20% polyethylene glycol (PEG)-6000. Unstressed seedlings without RA and PEG treatment were considered as control. At least three independent biological replicates of all four treatments namely, control, drought, RAinoculated and drought with RA inoculation were harvested at 72 h of stress treatment. Unstressed and stressed root samples were collected and snap frozen in liquid nitrogen and kept at −80 °C till further use.

2.5. Validation of miRNA expression using stem-loop quantitative real-time PCR (SL-qRT) The SL-qRT-PCR analysis was performed to validate the expression levels of selected conserved and novel miRNAs. Total RNA from all samples were reverse transcribed using TaqMan® MicroRNA Reverse Transcription kit (Applied Biosystems, Foster City, CA, USA) and realtime-PCR was done using DyNAmo Flash SYBR Green qPCR Kit (Thermo Scientific, USA) on a Stratagene Mx3000 P (Agilent, USA) Real-Time PCR system according to Jatan et al. (2018). SL-qRT-PCR primers were designed as described elsewhere (Chen et al., 2005; Supplementary Table S1). U6 snRNA was used as a reference gene for miRNAs expression data normalization (Jain et al., 2014). The relative expression of the miRNAs in different treated samples compared to the controls was calculated using 2−ΔΔCt method (Livak and Schmittgen, 2001).

2.2. Small RNA sequencing and data pre-processing Total RNA isolated from the control and treated root tissues of chickpea using Spectrum™ Plant Total RNA Kit (Sigma, USA) in three biological replicates each was pooled and used for library preparation using Illumina Truseq small RNA library preparation kit following the manufacturer's instructions (Illumina Technologies). The cluster generation and sequencing analysis of small RNA libraries was done by NxGenBio Life Sciences, Delhi, India using the Illumina Genome Analyzer platform (Illumina, USA). The small RNA sequence data (FASTQ files) generated has been deposited in the Sequence Read Archive (SRA) database under the accession number SRP137438. Raw reads produced by Illumina platform were pre-processed and read qualities were checked by using FastQC tool version 0.11.5. Adaptors and low-quality reads were removed by using adapter removal and filter tool of UEA sRNA workbench version 3.2 (Stocks et al., 2012).

2.6. Expression analysis of miRNAs targets by qRT-PCR To analyse the relative expression of respective target genes of all selected miRNAs qRT-PCR analysis was performed on Stratagene Mx3000 P (Agilent, USA) Real-Time PCR system. A chickpea glyceraldehydes 3-phosphate dehydrogenase (GAPDH) gene was used as an internal control for transcript normalization (Garg et al., 2010). The relative expression of the target genes was also determined using 2−ΔΔCt method (Livak and Schmittgen, 2001). Target gene-specific forward and reverse primers designing were done using IDT Primer Quest (https://eu.idtdna.com/PrimerQuest/Home/Index). List of target 218

Environmental and Experimental Botany 157 (2019) 217–227

R. Jatan et al.

miRNAs using miRCat. After processing, a total of 216 novel miRNAs were identified from all four libraries (Supplementary Table S5). The highest numbers of novel miRNAs were predicted from control (98) followed by drought (95), RA (82) and drought with RA (80) (Supplementary Table S6). However, 37, 36, 46 and 36 represented unique novel miRNAs from control, drought, RA, and drought with RA libraries, respectively. Out of 1139 miRNAs, 388 (34%) were existed in all four libraries and 391 (34%) were present in two or more libraries. A significant fraction (32%) of the miRNAs was unique and detected only from a particular library (Fig. 2A). The size distribution study of all miRNAs revealed the predominant (∼44%) presence of 21-nt-long miRNAs (Fig. 2B). Approximately, 26% represented 20-nt-long miRNAs and 25% miRNAs were of 18-19-nt in length. Remaining 5% miRNAs were of 22-24-nt in length. In addition, nucleotide sequence analysis of novel miRNAs showed that most predominant nucleotides was uridine at the 5′-end. Further, MFE of each predicted pre-miRNA secondary structure in chickpea was calculated to examine the accuracy of miRNA prediction. The average MFE of chickpea pre-miRNAs was found to be -58 kcal mol−1 and -42.54 kcal mol−1 for conserved and novel miRNAs, respectively. To support predicted miRNAs are true miRNAs, identification of complementary miRNA* sequences are important. Out of 216 novel miRNAs, the complementary miRNA* sequences for 29 novel miRNAs were detected.

gene-specific forward and reverse primers have been given in Supplementary Table S2. 2.7. Statistical analysis Relative expression data for miRNAs and their target genes from three independent biological replicates were calculated as the mean with standard error (mean ± SEM). Significant differences in variance between average values of control and treated plants, and comparison among averages was carried out using one-way analysis of variance (ANOVA) using Duncan's multiple range tests at p < .05 through Statistics Package for Social Science software version 16.0 (SPSS Inc./ IBM Corp. Chicago, USA). GraphPad Prism software (version 5.03, San Diego, California, USA) was used for illustrating the expression results. 3. Result 3.1. Small RNA sequencing To identify RA-responsive miRNAs from chickpea subjected to or not to drought stress, four sRNA libraries were generated from control, drought, RA and drought with RA-inoculated root tissues of chickpea. Illumina sequencing of all four libraries produced a total of 44,281,031 raw reads (Table 1). Approximately, 5, 11, 13, and 14 million raw reads were obtained from control, drought with RA, drought and RA-inoculated sRNAs libraries of chickpea root, respectively. After adapter trimming and selection of sequence length between 16 to 34 nucleotides, the total number of 1,068,810 reads from control, 2,449,400 reads from drought and 6,037,767 reads from RA and 1,773,927 reads from drought with RA library were obtained. After removal of tRNAs, rRNAs, snRNAs, and snoRNAs, the remaining 4,355,316 small RNA reads were utilised for miRNA identification. A total of 121,469 reads were mapped to chickpea genome. A high-quality sequencing data was indicated by the distribution of size of filtered sRNA sequence reads generated from all four libraries (Fig. 1A-D).

3.3. Identification of miRNA families All the identified conserved miRNAs were clustered into different miRNA families based on sequence similarity. Out of 923 conserved miRNAs, 914 were categorized into 24 families while rest nine miRNAs (miR363, miR384, miR5072, miR5083, miR5139, miR5168, miR5232, miR8155 and miR8175) did not show notable homology with other miRNAs (Fig. 3). Categorization of miRNAs families suggested that miR156, miR166 and miR167 were more frequently found miRNAs in chickpea genome. Seventeen miRNA families consisted of 11–80 miRNA members.

3.2. Identification and analysis of conserved and novel miRNAs Two approaches were used for the identification of chickpea miRNAs. In the first approach, the probable miRNA sequences that showed sequences homology with miRNAs of other plant species available in miRBase were considered as conserved miRNAs. A total of 923 conserved miRNAs were identified from all the libraries (Supplementary Table S3). Of 923 conserved miRNAs, 528, 758, 596, and 614 miRNAs were identified from control, drought, RA, and drought with RA (D + RA) libraries, respectively (Supplementary Table S4). Approximately 78% of conserved miRNAs were found in more than one library while 12, 80, 109 and 4 miRNAs were found to be unique in control, drought, RA, and drought with RA libraries, respectively. In the second approach, the putative small RNA reads from all four libraries (except those showed sequence homology with miRBase 21 entries) mapped on the genome of chickpea were used for prediction of novel

3.4. Distribution of novel miRNAs on chickpea chromosomes The genomic coordinates of novel miRNAs were used to determine the distribution of each novel miRNA on chickpea genome. Out of 216 novel miRNAs, the genomic coordinates for 171 miRNAs were retrieved which located on eight linkage groups, however, the remaining miRNAs were located on scaffold region of chickpea genome. Considerable clustering (tandem arrays) determined the localization of miRNAs on all eight chromosomes of chickpea (Fig. 4). The chromosomal distribution analysis of predicted novel miRNAs revealed that highest and lowest numbers of miRNAs were located on chromosome 4 and chromosome 8, respectively.

Table 1 Statistics of small RNA sequencing of all four libraries of chickpea. Properties

Control

Drought

RA-inoculated

Drought + RA-inoculated

Total number of raw reads Sequences remaining after adapter trimming Sequences remaining after length range filtering (16-35) filter low-complexity sequences filter invalid sequence filter by tRNA/rRNA Predicted conserved miRNAs Reads mapped to genome Predicted novel miRNAs

5404469 1806489 1068810 1068693 930454 587930 528 33346 98

13110618 3458037 2449400 2449214 2028584 1063175 758 32639 95

14542715 13180312 6037767 6037591 4893531 2316981 596 19878 82

11223229 2662393 1773927 1773879 1500651 1029818 614 35606 80

219

Environmental and Experimental Botany 157 (2019) 217–227

R. Jatan et al.

Fig. 1. (A–D). Size distribution of sRNA reads from chickpea roots. (A) Con (control), (B) D (drought), (C) RA (RA-inoculated), and (D) D + RA (drought with RAinoculated).

Fig. 2. (A–B). Predicted miRNAs in chickpea roots. (A) Number of identified miRNAs in all libraries (common), more than one library (multiple), and those unique to each library are given. (B) Nucleotide length distribution of miRNAs and the identification of the first nucleotide.

3.5. Targets of known and novel miRNAs To know the biological roles of miRNAs, it is essential to predict their target genes. A total of 4647 and 551 putative target genes for 762 conserved and 130 novel miRNAs, respectively were predicted from chickpea genome using psRNATarget server (Supplementary Tables S7 and S8). Out of 4647 and 551 predicted targets, 110 and 151 were identified as unique potential target genes for conserved and novel chickpea miRNAs, respectively.

3.6. GO enrichment analysis of miRNA target genes Fig. 3. Distribution of conserved miRNA families.

Blastx (version 2.2.29) search was performed against PlantRefSeq database for the predicted 5198 target genes of 1139 chickpea miRNAs and gene ontology (GO) enrichment analysis was also performed using 220

Environmental and Experimental Botany 157 (2019) 217–227

R. Jatan et al.

Fig. 4. Distribution of novel miRNAs on chickpea chromosomes.

Fig. 5. (A–C). Gene ontology (GO) enrichment analysis of the predicted targets of chickpea miRNAs. (A) biological processes, (B) cellular component and (C) molecular function.

221

Environmental and Experimental Botany 157 (2019) 217–227

R. Jatan et al.

Fig. 6. (A–D). The Venn diagram shows the numbers of common and unique differentially expressed miRNAs in chickpea. (A) and (B) conserved differentially expressed miRNAs, (C) and (D) novel differentially expressed miRNAs. (UP), up-regulated; (D), Down-regulated.

3.9. Comparative expression analysis of miRNAs and their target genes

UniProt database. GO analysis revealed that mostly the biological process terms of the target genes were involved in transcription (31%), regulation of transcription (30%), developmental process (9%), auxinactivated signaling pathways (7%) and cell division (4%) (Fig. 5A). The cellular component category included target genes which were involved in nucleus (80%) and mediator complex (9%) (Fig. 5B) while highly represented terms in molecular function of the target genes were involved in DNA binding (44%), metal ion binding (18%), transcription factor activity (7%), lipid binding (6%) and ATP binding (6%) (Fig. 5C).

In order to determine the correlation between expression patterns of miRNAs and their potential target genes, qRT-PCR analyses of predicted targets of above selected nine differentially expressed miRNAs was carried out under different treatments (Fig. 7A–I). Expression analysis was performed for mitogen-activated protein kinase kinase kinase 1like (target of miR156a), GAMYB-like (target of miR159a), auxin response factor 18-like (target of miR160a), homeobox-leucine zipper protein ATHB15-like (target of miR166 h), laccase-4-like (target of miR397a), basic blue protein-like (target of miR408), ethylene-responsive transcription factor 7-like (target of miR8175), uncharacterized protein with a putative protein kinase-like domain (target of Nov_miR71) and basic blue protein-like (target of Nov_miR104). Seven potential targets exhibited inverse correlations with the respective miRNAs expression patterns.

3.7. Differential expression analysis of miRNAs To know the putative roles of miRNAs, expression profiles of 923 conserved and 216 novel miRNAs from all four libraries was analysed. A total of 417, 431 and 332 conserved miRNAs were found to be differentially expressed under drought, RA, and drought with RA inoculation, respectively as compared to the control (Supplementary Tables S9) while 28, 59 and 35 novel miRNAs showed significant differential expression under drought, RA and drought with RA inoculation, respectively as compared to the control (Supplementary Tables S10). Out of 923 conserved miRNAs, 61, 149 and 83 miRNAs were significantly down-regulated in drought, RA and drought with RA, respectively in comparison to the control (Fig. 6A). However, 356, 282 and 249 miRNAs showed significant up-regulation in RA, drought and drought with RA-inoculated plants, respectively as compared to the control (Fig. 6B). Out of 216 novel miRNAs, 13, 35 and 12 miRNAs were significantly down-regulated in drought-treated, RA-inoculated, and drought with RA-inoculated roots when compared to the control, respectively (Fig. 6C). However, 15, 24 and 23 miRNAs were significantly up-regulated in drought, RA and drought with RA inoculation, respectively (Fig. 6D).

4. Discussion In recent years, high-throughput sequencing approaches have been widely utilised to identify miRNAs in numerous plants which play crucial roles in abiotic stress tolerance including response and tolerance to drought stress (Ding et al., 2013; Shuai et al., 2013; Zhang, 2015; Shriram et al., 2016). Nonetheless, PGPR-mediated drought stress tolerance through genome-wide modulation of miRNAs has yet not been explored. Therefore, in this study, we identified RA-responsive miRNAs in chickpea through high throughput sequencing that may regulate stress-responsive gene expression to withstand drought stress. Highest numbers of unique conserved and novel miRNAs were observed in RA-inoculated library as compared to the other libraries suggesting that the expression of miRNAs might be extensively regulated (either up-regulation or down-regulation) by RA-inoculation. RA has been previously reported to alter the expression of several TFs and other stress-responsive genes in Arabidopsis (Srivastva et al., 2012) and chickpea (Tiwari et al., 2016) as well as modulating the expression of miRNAs in chickpea upon drought and salt stress (Jatan et al., 2018). Taken together, RA may play important roles in the adaptation and restoration of cellular homeostasis in response to drought stress by the regulation of miRNAs expression. The difference in lengths of miRNAs (21–24 nt) produced is determined by the action of numerous DCL and AGO proteins (Cuperus et al., 2011; Jain et al., 2014). For example, 24nt long sRNAs produced by the action of DCL3 (Wu et al., 2011), and 24-nt miRNAs having an adenine at the 5′ end are characteristic features of AGO4 association (Mi et al., 2008). The 21-nt long miRNAs with 5′-uridine residue is a distinctive DCL1 cleavage and AGO1 association feature and has also been reported in most conserved miRNAs

3.8. Validation of expression profiles of selected conserved and novel miRNAs The expression patterns of seven conserved and two novel miRNAs were examined by using SL-qRT-PCR in all libraries (Fig. 7A–I). The results showed good concordance with those obtained by small RNA deep sequencing as evident from an overall correlation coefficient of 0.86 (Fig. 8). RA-inoculation resulted in highest up-regulation of miR408 (8.9 fold) and most markedly down-regulation of miR8175 (-3.0 fold) in comparison to the control. While Nov_miR104 was found to be highly up-regulated (5.3 fold) and Nov_miR71 showed most markedly down-regulation (-2.4 fold) as compared to the control. 222

Environmental and Experimental Botany 157 (2019) 217–227

R. Jatan et al.

Fig. 7. (A–I). Comparative expression analysis of conserved and novel miRNAs by RNASeq and SL-qRT-PCR along with their respective target genes by qRT-PCR in chickpea exposed to drought stress in the presence or absence of RA. Error bars represent the mean obtained from three biological replicates. Lowercase letters in the graph indicate significant differences (p-value ≤ 0.05) determined by Duncan’s multiple range test. Con, control; D, drought; RA, RA-inoculated; D + RA, drought with RA-inoculated.

report on chickpea (Jain et al., 2014). The predicted target genes of miR156a encoded a BTB/POZ (bric-a-brac, tramtrack, broad complex/ Pox virus and zinc finger) domain and TAZ (Transcriptional adapter zinc binding) domain-containing protein, mitogen-activated protein kinase kinase kinase 1-like (MAPKKK-1 like) and Squamosa PromoterBinding-Like (SPL) protein. The BTB /POZ domain is a multifaceted protein-protein interaction motif that contributes towards a wide range of cellular functions, including transcriptional regulation, and protein ubiquitination in plants (Stogios et al. 2005). Similarly, SPL family of transcription factors (TFs) is also involved in the regulation of fundamental aspects of growth and development of plants such as leaf initiation rate, branching, vegetative phase change and flowering time (Preston and Hileman, 2013). While, MAPKKK-1 like has been reported to be involved in abiotic and biotic stress signalling pathways and stress response in rice (Raghuram et al., 2014). Taken together, the most frequent occurrence of miR156 members in chickpea genome indicates their versatile role and involvement in complex regulatory mechanisms

(Mi et al., 2008; Breakfield et al., 2012). In this study, the size distribution analysis of filtered sRNA sequence reads revealed that approximately 53% of sRNA sequence reads were within the range of 21–24 nt (Fig. 1), as expected, which is a characteristic feature of the DCL cleavage products reported in earlier studies (Breakfield et al., 2012; Hwang et al., 2013) including a study on chickpea (Jain et al., 2014). MFE usually determines the stability of RNA secondary structure and has been utilised for miRNA prediction (Llave et al., 2002; Zhang et al., 2009). A lower MFE value indicates highly stable RNA secondary structure, which is a characteristic feature of miRNAs (Bonnet et al., 2004). The average MFE of chickpea pre-miRNAs was found to be lower (-58 kcal mol−1 for conserved and -42.54 kcal mol−1 for novel miRNAs. Similar findings have been reported in previous studies too (Thakur et al., 2011; Jain et al., 2014). In our study, miR156 and miR166 constituted the largest miRNA families in chickpea with 131 and 121 members, respectively. The findings are in accordance to a previous 223

Environmental and Experimental Botany 157 (2019) 217–227

R. Jatan et al.

by both miR397 and miR408. Apart from conserved miRNA families, various novel miRNAs were also predicted to target Ras-related protein RABB1c, aspartate-tRNA ligase, serine/threonine-protein kinase, calcium-transporting ATPase, peroxidase superfamily protein, ATP sulfurylase, potassium channel KAT1, Nodulation-signaling pathway 2 protein (NSP2) and numerous TFs protein-encoding genes. The previous study reported that miR171c and miR397 targets NSP2 and laccase genes respectively, involved in root nodulation and maintenance of copper homeostasis in Lotus japonicus root nodules (De Luis et al., 2012). However, in our study, NSP2 protein-coding gene was also found to be targeted by a novel miRNA36 (Nov_miR36) and a laccase coding gene was predicted to be targeted by Nov_miR79 and Nov_miR104 family members. The presence of Nov_miR36 and Nov_miR104 specifically in RA-inoculated library suggested that RA might be involved in root nodulation and copper homeostasis in chickpea as well. ATP sulfurulases involved in the assimilation of sulphate and predicted to be a target of Nov_miR115 was found to be differentially expressed in both drought with RA and RA-inoculated libraries indicating the possible involvement of RA in sulphate assimilation. Further, functional categorization of the proteins encoded by predicted target genes of miRNAs demonstrated their involvement in biological, cellular and molecular functions. The highest GO terms were observed for transcription, regulation of transcription in the biological process, terms for the nucleus in the cellular component, and terms for DNA binding in molecular function categories indicating the reliability of our data as miRNAs are basically involved in regulation of target genes in nucleus. In addition, this study also throws insights into the role of RA in growth promotion of chickpea under drought stress. As for example, a ubiquitin-conjugating E2 enzyme (UBC24) encoded by PHOSPHATE2 (PHO2) involved in protein degradation, is a target of miR399 which down-regulates its expression under Pi-deficient condition, consequently stimulating Pi uptake and translocation of Pi from root to shoot to maintain Pi homeostasis (Aung et al., 2006; Bari et al., 2006). In our study, miR399 was found to be upregulated only in drought and drought with RA-inoculated conditions and not in RA-inoculated plants suggesting that RA may be involved in maintaining phosphate homeostasis in chickpea through its phosphate solubilisation activity. In our study, a new miRNA, miR5232 was predicted to target calcium-transporting ATPase which is reported to be involved in growth, development and adaptation to adverse environmental conditions (Frey et al., 2012; Huda et al., 2013). miR5232 was found to be differentially expressed only in drought and drought with RA-inoculated conditions suggesting the crucial role of RA in growth, development and drought stress alleviation in chickpea. Further, quantitative real-time PCR analysis of selected nine miRNAs and their respective targets revealed differential expression patterns under all the treatments studied. Moreover, the SL-qRT-PCR expression patterns of selected miRNAs under different treatments were largely same to that obtained from Illumina sequencing suggesting that the drought and RA-responsive miRNAs identified in this study are reliable. A GAMYB-like TF has been predicted to be targeted by miR159a and has been shown to play crucial roles in hormone signaling, stomatal regulation and drought stress response (Baldoni et al., 2015). In our study, miR159 and GAMYB-like TF showed a negative correlation under all treatments as compared to the control indicating that this gene might play a very crucial role in drought stress response/ tolerance in chickpea. Auxin-responsive factor (ARF) genes namely, ARF10, ARF16, ARF17,) and ARF18 have been reported to be the targets of miR160 which leadto stress adaptation by the regulation of numerous physiological and molecular processes in various plants (Wang et al., 2005; Guilfoyle and Hagen, 2007; Candar-Cakir et al., 2016; Sun et al., 2017). In this study, miR160a was found to be significantly up-regulated in RA-inoculated plants and an inverse correlation in the expression of miR160a and ARF18 indicated the involvement of RA in ameliorating drought stress by regulating miR160a and its target.

Fig. 8. Correlation of gene expression results obtained from small RNA-seq and SLqRT-PCR analysis for nine selected miRNAs from four libraries. Each data point displays the log2 normalized expression level obtained from small RNAseq (y axis) and SL-qRT-PCR (x axis) analyses.

that might play important roles in growth and development, transcriptional regulation and response to drought stress in chickpea upon RA-inoculation. The physical mapping of novel miRNA families on chickpea genome suggested the existence of more than one group of miRNAs of the same family on one or more chromosomes (Fig. 4). Such type of tandem clustering of numerous miRNAs families has been previously reported from different plants including Arabidopsis, rice and Medicago (JonesRhoades and Bartel, 2004; Cui et al., 2009; Lelandais-Briere et al., 2009). Based on sequence homology, we assumed that tandem and/or segmental duplication lead to the expansion of miRNA gene families in chickpea to regulate the expression of genes in a remarkably dosagedependent and spatiotemporal manner as reported in several plants earlier (Li and Mao, 2007; Cuperus et al., 2011). Majority of the predicted target genes of various conserved and novel miRNAs were found to encode TFs such as SPB, NAC [No apical meristem (NAM); Arabidopsis transcription activation factor (ATAF); Cup-shaped cotyledon (CUC)], GRF (Growth regulating factor), Zinc finger and MYB (Myeloblastosis), class III HD-ZIP (Class III homeodomain-leucine zipper), auxin responsive factor (ARFs), AP2-EREBP [AP2 (APETALA2); EREBPs (ethylene-responsive element binding proteins)], GRAS [Gibberellic-acid insensitive (GAI), Repressor of GAI (RGA) and Scarecrow (SCR)], and NF-Y (Nuclear transcription factor Y) TF families in accordance to previous studies (Bartel, 2004; Jeong and Zhai, 2011; Breakfield et al., 2012; Hwang et al., 2013; Jatan et al., 2018). In addition to the above mentioned targets, novel miRNAs were also predicted to target various other protein-coding genes including kinases and transporter proteins. Interestingly, numerous target genes of conserved (31%) and novel miRNAs (42%) predicted in our study were genes that encoded uncharacterised proteins and genes with unknown functions indicating the existence of novel and complex posttranscriptional gene regulatory networks that may be operating in chickpea upon RA-inoculation under drought stress. The analysis of target inhibition process in this study revealed that more than 85% targets for conserved and novel miRNAs were controlled through cleavage and the remaining by translational inhibition as previously reported in plants including chickpea (Mallory and Vaucheret, 2010; Rogers and Chen, 2013; Jain et al., 2014). We have found that almost all identified miRNAs can regulate more than one transcript, however some transcripts can be targeted by more than one miRNA also which is in accordance to a previous study in Dimocarpus longan (Lin and Lai, 2013). As for example, miR156 was predicted to target three genes described earlier, miR393 targeted both Auxin Signaling F-Box 2 and Transport Inhibitor Response 1, and miR408 targeted basic blue protein, laccase and uclacyanin. While miR159 and miR319 were predicted to target MYB81, GAMYB and TCP TFs, and laccase gene was targeted 224

Environmental and Experimental Botany 157 (2019) 217–227

R. Jatan et al.

Fig. 9. Schematic representation of the role of miRNAs in response to drought stress in chickpea in the presence of Pseudomonas putida RA.

plantacyanin expression and an inducer of DREB (dehydration-responsive element binding) gene under drought stress (Hajyzadeh et al., 2015). In this study, miR408b and Nov_miR104 were found to be significantly up-regulated and BBP-like gene significantly down-regulated in RA-inoculated plants as compared to the control indicating that RAinoculation affects the accumulation of copper by regulating miR408b, Nov_miR104 and BBP-like genes, thus maintaining elemental balance and redox homeostasis. Ethylene-Responsive Factor (ERF) TF belongs to the APETALA2/ ERF superfamily that regulates plant growth and development as well as biotic and abiotic stress tolerance (Lata et al., 2011; Lata and Prasad, 2011; Mehrnia et al., 2013). ERF7 has been reported to play important roles in response to drought stress in sugarcane (Vargas et al., 2014). In our study, miR8175 was found to be significantly down-regulated and inversely regulate the expression of ERF7 under all treatments indicating the involvement of miR8175-ERF7 in drought stress response and tolerance in chickpea. The predicted target gene of Nov_miR71 encodes an uncharacterized protein with a putative protein kinase-like (PK-like) domain. Protein kinases have been reported to be involved in ionic and osmotic homeostasis signaling pathways as well as drought stress response in plants (Deng et al., 2017). In our study, an opposite expression pattern of Nov_miR71 with its target gene was found in drought with RA-inoculated plants. However, the functions of uncharacterized protein with a putative PK-like domain remains elusive and need to be explored further to understand its role in PGPR-mediated stress response and tolerance in chickpea.

miR166 is an important drought responsive miRNA that targets the HD-ZIP III TF family members involved in vascular development and stress tolerance (Kim, 2005; Ding et al., 2009). In our study, miR166 h was found to be significantly up-regulated under drought stress which is in accordance to previous reports in barley, Medicago and Saccharum under drought stress (Trindade et al., 2010; Kantar et al., 2011; Gentile et al., 2015). In RA-inoculated plants, miR166h and ATHB15 showed opposite expression pattern suggesting the role of RA in alleviating drought stress and promoting growth in inoculated plants. In our study, RA-inoculated plants showed down-regulation of miR156a and an upregulation of MAPKKK-1 like gene, indicating that RA might be involved in growth and development of chickpea by the regulation of miR156a-MAPKKK module. The same is in accordance with genomewide expression profiles of two banana genotypes that highlighted the involvement of MAPKK and MAPKKK genes in banana tissue and fruit development, fruit ripening and responses to abiotic stresses (Wang et al., 2017). The overexpression of miR397a has been reported to down-regulate the expression of laccase genes leading to reduced lignin deposition in both transgenic P. trichocarpa (Lu et al., 2013), and transgenic Arabidopsis (Wang et al., 2014). In Arabidopsis, overexpression of miR397a also increased number and length of siliques resulting in higher seed numbers demonstrating regulation of both lignin biosynthesis and grain yield by this miRNA(Wang et al., 2014). Accordingly, in this study the expression of laccase-4-like gene showed significant down-regulation in drought, RA-inoculated and drought with Ra-inoculated plants while miR397a was found to be either basal or significantly up-regulated at all conditions analyzed suggesting RA-mediated miR397-laccase gene regulation might be crucial for improved seed yield and reduced lignin content in chickpea too. miR408 family members have been reported to regulate the expression of plantacyanin-like (basic blue) proteins in Arabidopsis and wheat (Sunkar and Zhu, 2004; Yao et al., 2007). In our study, other than miR408a, a novel miRNA104 (Nov_miR104) was also found to target a basic blue protein-like (BBP-like) gene, although there exists a single nucleotide difference between these two miRNAs. The overexpression of miR408 has been reported to affect the accumulation of copper in chickpea and was established as a negative regulator of

5. Conclusion Overall, we are for the first time reporting on the genome-wide profiling of miRNAs of chickpea after colonization by P. putida RA under drought stress using small RNA sequencing approach. Though drought tolerance is a complex phenomenon, RA-mediated regulation of expression of miRNAs and their target genes helped us to speculate that the adaptive responses associated to normal physiology, growth and development of plants under drought stress are modulated by RAinoculation. A working hypothesis in chickpea for the mechanism of 225

Environmental and Experimental Botany 157 (2019) 217–227

R. Jatan et al.

miRNA-mediated drought tolerance under the influence of RA has also been elaborated (Fig. 9). Negative expression pattern of miRNAs and their target genes indicated the involvement in the stress regulatory pathways, which control/regulate drought stress response during drought stress in chickpea. This study thus revealed that the RA-inoculation might play crucial roles in the modulation of miRNAs and their targets in drought stress amelioration and their utilization in crop improvement programs.

Stougaard, J., Voinnet, O., 2012. Two microRNAs linked to nodule infection and nitrogen-fixing ability in the legume Lotus japonicus. Plant Physiol. 160, 2137–2154. Deng, M., Cao, Y., Zhao, Z., Yang, L., Zhang, Y., Dong, Y., Fan, G., 2017. Discovery of MicroRNAs and their target genes related to drought in Paulownia "yuza 1" by Highthroughput sequencing. Int. J. Genomics 3674682. Ding, D., Zhang, L., Wang, H., Liu, Z., Zhang, Z., Zheng, Y., 2009. Differential expression of miRNAs in response to salt stress in maize roots. Ann. Bot. 103, 29–38. Ding, Y., Tao, Y., Zhu, C., 2013. Emerging roles of microRNAs in the mediation of drought stress response in plants. J. Exp. Bot. 64, 3077–3086. Ferdous, J., Sanchez-Ferrero, J.C., Langridge, P., Milne, L., Chowdhury, J., Brien, C., Tricker, P.J., 2017. Differential expression of microRNAs and potential targets under drought stress in barley. Plant Cell Environ. 40, 11–24. Frey, N.F.D., Mbengue, M., Kwaaitaal, M., Nitsch, L., Altenbach, D., Haweker, H., Lonzano-Duran, R., Njo, M.F., Beeckman, T., Huettel, B., Borst, J.W., 2012. Plasma membrane calcium ATPases are important components of receptor-mediated signalling in plant immune responses and development. Plant Physiol. 111. Garg, R., Sahoo, A., Tyagi, A.K., Jain, M., 2010. Validation of internal control genes for quantitative gene expression studies in chickpea (Cicer arietinum L.). Biochem. Biophys. Res. Commun. 396, 283–288. Gentile, A., Dias, L.I., Mattos, R.S., Ferreira, T.H., Menossi, M., 2015. MicroRNAs and drought responses in sugarcane. Front. Plant Sci. 6, 58. Griffiths-Jones, S., Saini, H.K., van Dongen, S., Enright, A.J., 2008. miRBase: tools for microRNA genomics. Nucleic Acids Res. 36, D154–D158. Guilfoyle, T.J., Hagen, G., 2007. Auxin response factors. Curr. Opin. Plant Biol. 10, 453–460. Hajyzadeh, M., Turktas, M., Khawar, K.M., Unver, T., 2015. miR408 overexpression causes increased drought tolerance in chickpea. Gene 555, 186–193. Huda, K.M.K., Banu, M.S.A., Tuteja, R., Tuteja, N., 2013. Global calcium transducer Ptype Ca2+-ATPases open new avenues for agriculture by regulating stress signalling. J. Exp. Bot. 64 (11), 3099–3109. Hwang, D.G., Park, J.H., Lim, J.Y., et al., 2013. The hot pepper (capsicum annuum) microRNA transcriptome reveals novel and conserved targets: a foundation for understanding MicroRNA functional roles in hot pepper. PLoS One 8, e64238. Jain, M., Misra, G., Patel, R.K., Priya, P., Jhanwar, S., Khan, A.W., Shah, N., Singh, V.K., Garg, R., Jeena, G., Yadav, M., 2013. A draft genome sequence of the pulse crop chickpea (Cicer arietinum L.). The Plant J. 1;74 (5), 715–729. Jain, M., Chevala, V.V., Garg, R., 2014. Genome-wide discovery and differential regulation of conserved and novel microRNAs in chickpea via deep sequencing. J. Exp. Bot. 65, 5945–5958. Jatan, R., Chauhan, P.S., Lata, C., 2018. Pseudomonas putida modulates the expression of miRNAs and their target genes in response to drought and salt stresses in chickpea (Cicer arietinum L.). Genomics. https://doi.org/10.1016/j.ygeno.2018.01.007. Jeong, D.H., Zhai, Park S., 2011. Massive analysis of rice small RNAs: mechanistic implications of regulated microRNAs and variants for differential target RNA cleavage. The Plant Cell 23, 4185–4207. Jones-Rhoades, M.W., Bartel, D.P., 2004. Computational identification of plant microRNAs and their targets, including a stress-induced miRNA. Mol. Cell 14, 787–799. Kantar, M., Lucas, S.J., Budak, H., 2011. miRNA expression patterns of triticum dicoccoides in response to shock drought stress. Planta 233, 471–484. Khandal, H., Parween, S., Roy, R., Meena, M.K., Chattopadhyay, D., 2017. MicroRNA profiling provides insights into post-transcriptional regulation of gene expression in chickpea root apex under salinity and water deficiency. Scientific Rep. 7, 4632. Khodadadi, M., 2013. Effect of drought stress on yield and water relative content in chickpea. Int. J. Agron. Plant Product. 6, 1168–1172. Kim, V.N., 2005. MicroRNA biogenesis: coordinated cropping and dicing. Nat. Rev. Mol. Cell Biol. 6, 376–385. Kohli, D., Joshi, G., Deokar, A.A., Bhardwaj, A.R., Agarwal, M., Katiyar-Agarwal, S., Srinivasan, R., Jain, P.K., 2014. Identification and characterization of Wilt and salt stress-responsive microRNAs in chickpea through high-throughput sequencing. PLoS One 9, e108851. Lata, C., Prasad, M., 2011. Role of DREBs in regulation of abiotic stress responses in plants. J. Exp. Bot. 62, 4731–4748. Lata, C., Bhutty, S., Bahadur, R.P., Majee, M., Prasad, M., 2011. Association of an SNP in a novel DREB2-like gene SiDREB2 with stress tolerance in foxtail millet [setaria italica (L.)]. J. Exp. Bot. 62, 3387–3401. Lata, C., Muthamilarasan, M., Prasad, M., 2015. Drought Stress Responses and Signal Transduction in Plants. In: Elucidation of Abiotic Stress Signaling in Plants. Springer, pp. 195–225. Lelandais-Briere, C., Naya, L., Sallet, E., Calenge, F., Frugier, F., Hartmann, C., Gouzy, J., Crespi, M., 2009. Genome-wide Medicago truncatula small RNA analysis revealed novel microRNAs and isoforms differentially regulated in roots and nodules. The Plant Cell 21, 2780–2796. Li, H., Durbin, R., 2009. Fast and accurate short read alignment with burrows-wheeler transform. Bioinformatics 25, 1754–1760. Li, A., Mao, L., 2007. Evolution of plant microRNA gene families. Cell Res. 17, 212–218. Li, B., Duan, H., Li, J., Deng, X.W., Yin, W., Xia, X., 2013. Global identification of miRNAs and targets in Populus euphratica under salt stress. Plant Mol. Biol. 81, 525–539. Lin, Y., Lai, Z., 2013. Comparative analysis reveals dynamic changes in miRNAs and their targets and expression during somatic embryogenesis in longan (Dimocarpus longan Lour.). PLoS One 8, e60337. Liu, F., Wang, W., Sun, X., Liang, Z., Wang, F., 2015. Conserved and novel heat stressresponsive microRNAs were identified by deep sequencing in Saccharina japonica (Laminariales, Phaeophyta). Plant Cell Environ. 38, 1357–1367. Liu, M., Yu, H., Zhao, G., Huang, Q., Lu, Y., Ouyang, B., 2017. Profiling of droughtresponsive microRNA and mRNA in tomato using high-throughput sequencing. BMC

Author contribution statement C.L. conceived and designed the experiment. R.J. and S.T. conducted the experiments. R.J., M.H.A. and C.L. analyzed the data. R.J. and C.L. wrote the manuscript. All authors read and approved the manuscript. Conflict of interest The authors declare that they have no conflict of interest. Acknowledgments CL acknowledges SERB Women Excellence Award [Grant No. SB/ WEA-05/2014] from Science & Engineering Research Board (SERB), Government of India, New Delhi. RJ acknowledges Junior Research Fellowship (Fellow No. DBT/2015/NBRI/322) from Department of Biotechnology, Government of India. Authors are highly thankful to Dr. Puneet S. Chauhan for providing the bacterial strain. Authors are hig Appendix A. Supplementary data Supplementary material related to this article can be found, in the online version, at doi:https://doi.org/10.1016/j.envexpbot.2018.10. 003. References Ahuja, I., de Vos, R.C., Bones, A.M., Hall, R.D., 2010. Plant molecular stress responses face climate change. Trends Plant Sci. 15, 664–674. Aravind, J., Rinku, S., Pooja, B., Shikha, M., Kaliyugam, S., Mallikarjuna, M.G., Kumar, A., Rao, A.R., Nepolean, T., 2017. Identification, characterization, and functional validation of drought-responsive MicroRNAs in subtropical maize inbreds. Front. Plant Sci. 8, 941. Aung, K., Lin, S.I., Wu, C.C., Huang, Y.T., Su, C.L., Chiou, T.J., 2006. pho2, a phosphate overaccumulator, is caused by a nonsense mutation in a microRNA399 target gene. Plant physiol 141, 1000–1011. Baldoni, E., Genga, A., Cominelli, E., 2015. Plant MYB transcription factors: their role in drought response mechanisms. Int. J. Mol. Sci. 16, 15811–15851. Barakat, A., Sriram, A., Park, J., Zhebentyayeva, T., Main, D., Abbott, A., 2012. Genome wide identification of chilling responsive microRNAs in Prunus persica. BMC Genomics 13, 481. Bari, R., Pant, B.D., Stitt, M., Scheible, W.R., 2006. PHO2, microRNA399, and PHR1 define a phosphate-signaling pathway in plants. Plant Physiol. 141, 988–999. Bartel, D.P., 2004. MicroRNAs: genomics, biogenesis, mechanism, and function. Cell 116, 281–297. Bonnet, E., Wuyts, J., Rouze, P., Van de Peer, Y., 2004. Evidence that microRNA precursors, unlike other non-coding RNAs, have lower folding free energies than random sequences. Bioinformatics 20, 2911–2917. Breakfield, N.W., Corcoran, D.L., Petricka, J.J., Shen, J., Sae-Seaw, J., Rubio-Somoza, I., Weigel, D., Ohler, U., Benfey, P.N., 2012. High-resolution experimental and computational profiling of tissue-specific known and novel miRNAs in Arabidopsis. Genome Res. 22, 163–176. Candar-Cakir, B., Arican, E., Zhang, B., 2016. Small RNA and degradome deep sequencing reveals drought-and tissue-specific micrornas and their important roles in droughtsensitive and drought-tolerant tomato genotypes. Plant Biotechnol. J. 14, 1727–1746. Chen, C., Ridzon, D.A., Broomer, A.J., et al., 2005. Real-time quantification of microRNAs by stem-loop RT-PCR. Nucleic Acids Res. 33, e179. Cui, X., Xu, S.M., Mu, D.S., Yang, Z.M., 2009. Genomic analysis of rice microRNA promoters and clusters. Gene 431, 61–66. Cuperus, J.T., Fahlgren, N., Carrington, J.C., 2011. Evolution and functional diversification of MIRNA genes. The Plant Cell 23, 431–442. Dai, X., Zhao, P.X., 2011. psRNATarget: a plant small RNA target analysis server. Nucleic Acids Res. 39, W155–159. De Luis, A., Markmann, K., Cognat, V., Holt, D.B., Charpentier, M., Parniske, M.,

226

Environmental and Experimental Botany 157 (2019) 217–227

R. Jatan et al.

Sunkar, R., Zhu, J.K., 2004. Novel and stress-regulated microRNAs and other small RNAs from Arabidopsis. The Plant Cell 16, 2001–2019. Talebi, R., Ensafi, M.H., Baghebani, N., Karami, E., Mohammadi, K., 2013. Physiological responses of chickpea (Cicer arietinum) genotypes to drought stress. Environ. Exp. Biol. 11, 9–15. Tiwari, S., Lata, C., Chauhan, P.S., Nautiyal, C.S., 2016. Pseudomonas putida attunes morphophysiological, biochemical and molecular responses in Cicer arietinum L. during drought stress and recovery. Plant Physiol. Biochem. 99, 108–117. Tiwari, S., Lata, C., Chauhan, P.S., Prasad, V., Prasad, M., 2017. A functional genomic perspective on drought signalling and its crosstalk with phytohormone-mediated signalling pathways in plants. Curr. Genomics 18 (6), 469–482. Trindade, I., Capitao, C., Dalmay, T., Fevereiro, M.P., Santos, D.M., 2010. miR398 and miR408 are up-regulated in response to water deficit in Medicago truncatula. Planta 231, 705–716. Vargas, L., Santa Brigida, A.B., Mota Filho, J.P., et al., 2014. Drought tolerance conferred to sugarcane by association with Gluconacetobacter diazotrophicus: a transcriptomic view of hormone pathways. PLoS One 9, e114744. Varshney, R.K., Song, C., Saxena, R.K., et al., 2013. Draft genome sequence of chickpea (Cicer arietinum) provides a resource for trait improvement. Nat. Biotechnol. 31, 240–246. Voinnet, O., 2009. Origin, biogenesis, and activity of plant microRNAs. Cell 136, 669–687. Voorrips, R.E., 2002. MapChart: software for the graphical presentation of linkage maps and QTLs. J. Heredity 93, 77–78. Wang, J.W., Wang, L.J., Mao, Y.B., Cai, W.J., Xue, H.W., Chen, X.Y., 2005. Control of root cap formation by MicroRNA-targeted auxin response factors in Arabidopsis. The Plant Cell 17, 2204–2216. Wang, C.Y., Zhang, S., Yu, Y., et al., 2014. MiR397b regulates both lignin content and seed number in Arabidopsis via modulating a laccase involved in lignin biosynthesis. Plant Biotechnol. J. 12, 1132–1142. Wang, L., Hu, W., Tie, W., et al., 2017. The MAPKKK and MAPKK gene families in banana: identification, phylogeny and expression during development, ripening and abiotic stress. Sci. Rep. 7, 1159. Wu, L., Zhou, H., Zhang, Q., Zhang, J., Ni, F., Liu, C., Qi, Y., 2011. DNA methylation mediated by a microRNA pathway. Mol. Cell 38, 465–475. Wu, J., Wang, L., Wang, S., 2017. MicroRNAs associated with drought response in the pulse crop common bean (Phaseolus vulgaris L.). Gene 628, 78–86. Xia, H., Zhang, L., Wu, G., Fu, C., Long, Y., Xiang, J., Gan, J., Zhou, Y., Yu, L., Li, M., 2016. Genome-Wide identification and characterization of MicroRNAs and target genes in Lonicera japonica. PLoS One 11, e0164140. Xue, M., Yi, H., Wang, H., 2018. Identification of miRNAs Involved in SO 2 Preservation in Vitis vinifera L. by Deep Sequencing. Environmental and Experimental Botany. Yang, X., Liu, F., Zhang, Y., Wang, L., Cheng, Y.F., 2017. Cold-responsive miRNAs and their target genes in the wild eggplant species Solanum aculeatissimum. BMC Genomics 18, 1000. Yao, Y., Guo, G., Ni, Z., Sunkar, R., Du, J., Zhu, J.K., Sun, Q., 2007. Cloning and characterization of microRNAs from wheat (Triticum aestivum L.). Genome Biol. 8, R96. Yin, Z., Han, X., Li, Y., Wang, J., Wang, D., Wang, S., Fu, X., Ye, W., 2017. Comparative analysis of Cotton small RNAs and their target genes in response to Salt stress. Genes (Basel) 8. Zhang, B., 2015. MicroRNA: a new target for improving plant tolerance to abiotic stress. J. Exp. Bot. 66, 1749–1761. Zhang, J., Xu, Y., Huan, Q., Chong, K., 2009. Deep sequencing of Brachypodium small RNAs at the global genome level identifies microRNAs involved in cold stress response. BMC Genomics 10, 449.

Genomics 18, 481. Livak, K.J., Schmittgen, T.D., 2001. Analysis of relative gene expression data using realtime quantitative PCR and the 2(-Delta Delta C(T)) method. Methods 25, 402–408. Llave, C., Xie, Z., Kasschau, K.D., Carrington, J.C., 2002. Cleavage of scarecrow-like mRNA targets directed by a class of Arabidopsis miRNA. Science 297, 2053–2056. Lu, S., Li, Q., Wei, H., et al., 2013. Ptr-miR397a is a negative regulator of laccase genes affecting lignin content in Populus trichocarpa. Proc. Natl Acad. Sci. U. S. A. 110, 10848–10853. Ma, X., Xin, Z., Wang, Z., Yang, Q., Guo, S., Guo, X., Cao, L., Lin, T., 2015. Identification and comparative analysis of differentially expressed miRNAs in leaves of two wheat (Triticum aestivum L.) genotypes during dehydration stress. BMC Plant Biol. 15, 21. Macovei, A., Tuteja, N., 2012. microRNAs targeting DEAD-box helicases are involved in salinity stress response in rice (Oryza sativa L.). BMC Plant Biol. 12, 183. Mallory, A., Vaucheret, H., 2010. Form, function, and regulation of ARGONAUTE proteins. The Plant Cell 110 pp.tpc. Mehrnia, M., Balazadeh, S., Zanor, M.I., Mueller-Roeber, B., 2013. EBE, an AP2/ERF transcription factor highly expressed in proliferating cells, affects shoot architecture in Arabidopsis. Plant Physiol. 162, 842–857. Meyers, B.C., Axtell, M.J., Bartel, B., et al., 2008. Criteria for annotation of plant microRNAs. The Plant Cell 20, 3186–3190. Mi, S., Cai, T., Hu, Y., et al., 2008. Sorting of small RNAs into Arabidopsis argonaute complexes is directed by the 5’ terminal nucleotide. Cell 133, 116–127. Preston, J.C., Hileman, L.C., 2013. Functional evolution in the plant SQUAMOSAPROMOTER BINDING PROTEIN-LIKE (SPL) gene family. Front. Plant Sci. 4, 80. Raghuram, B., Sheikh, A.H., Sinha, A.K., 2014. Regulation of MAP kinase signaling cascade by microRNAs in Oryza sativa. Plant Signal. Behav. 9, e972130. Robinson, M.D., McCarthy, D.J., Smyth, G.K., 2010. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139–140. Rogers, K., Chen, X., 2013. Biogenesis, turnover, and mode of action of plant microRNAs. The Plant Cell 25, 2383–2399. Shivhare, R., Lata, C., 2017. Exploration of genetic and genomic resources for abiotic and biotic stress tolerance in Pearl Millet. Front. Plant Sci. 7, 2069. Shriram, V., Kumar, V., Devarumath, R.M., Khare, T.S., Wani, S.H., 2016. MicroRNAs as potential targets for abiotic stress tolerance in plants. Front. Plant Sci. 7, 817. Shuai, P., Liang, D., Zhang, Z., Yin, W., Xia, X., 2013. Identification of drought-responsive and novel Populus trichocarpa microRNAs by high-throughput sequencing and their targets using degradome analysis. BMC Genomics 14, 233. Srivastava, S., Chaudhry, V., Mishra, A., Chauhan, P.S., Rehman, A., Yadav, A., Tuteja, N., Nautiyal, C.S., 2012. Gene expression profiling through microarray analysis in Arabidopsis thaliana colonized by Pseudomonas putida MTCC5279, a plant growth promoting rhizobacterium. Plant Signal. Behav. 7, 235–245. Srivastava, S., Zheng, Y., Kudapa, H., Jagadeeswaran, G., Hivrale, V., Varshney, R.K., Sunkar, R., 2015. High throughput sequencing of small RNA component of leaves and inflorescence revealed conserved and novel miRNAs as well as phasiRNA loci in chickpea. Plant Sci. 235, 46–57. Stocks, M.B., Moxon, S., Mapleson, D., Woolfenden, H.C., Mohorianu, I., Folkes, L., Schwach, F., Dalmay, T., Moulton, V., 2012. The UEA sRNA workbench: a suite of tools for analysing and visualizing next generation sequencing microRNA and small RNA datasets. Bioinformatics 28, 2059–2061. Sun, R., Li, C., Zhang, J., Li, F., Ma, L., Tan, Y., Wang, Q., Zhang, B., 2017. Differential expression of microRNAs during fiber development between fuzzless-lintless mutant and its wild-type allotetraploid cotton. Sci. Rep. 7, 3. Sunkar, R., Li, Y.F., Jagadeeswaran, G., 2012. Functions of microRNAs in plant stress responses. Trends Plant Sci. 17, 196–203.

227