Integrative small RNA and transcriptome analysis provides insight into key role of miR408 towards drought tolerance response in cowpea

Drought stress response studies and overexpression of vun-miR408 proved it to be essential for abiotic stress tolerance in cowpea. Small RNA and transcriptome sequencing of an elite high-yielding drought-tolerant Indian cowpea cultivar, Pusa Komal revealed a differential expression of 198 highly conserved, 21 legume-specific, 14 less-conserved, and 10 novel drought-responsive microRNAs (miRNAs) along with 3391 (up-regulated) and 3799 (down-regulated) genes, respectively, in the leaf and root libraries. Among the differentially expressed miRNAs, vun-miR408-3p, showed an up-regulation of 3.53-log2-fold change under drought stress. Furthermore, laccase 12 (LAC 12) was identified as the potential target of vun-miR408-3p using 5′ RNA ligase-mediated rapid amplification of cDNA ends. The stable transgenic cowpea lines overexpressing artificial vun-miR408-3p (OX-amiR408) displayed enhanced drought and salinity tolerance as compared to the wild-type plants. An average increase of 30.17% in chlorophyll, 26.57% in proline, and 27.62% in relative water content along with lesser cellular H2O2 level was observed in the transgenic lines in comparison with the wild-type plants under drought stress. Additionally, the scanning electron microscopic study revealed a decrease in the stomatal aperture and an increase in the trichome density in the transgenic lines. The expression levels of laccase 3 and laccase 12, the potential targets of miR408, related to lipid catabolic processes showed a significant reduction in the wild-type plants under drought stress and the transgenic lines, indicating the regulation of lignin content as a plausibly essential trait related to the drought tolerance in cowpea. Taken together, this study primarily focused on identification of drought-responsive miRNAs and genes in cowpea, and functional validation of role of miR408 towards drought stress response in cowpea.


Introduction
Drought stress is considered as the most destructive and an inescapable reality gradually prevalent in most parts of the world further intensified with the unavoidable social and ecological disturbances, leading to higher occurrence of agriculturally unfavoured parched land areas, reduction in vegetation, increased wild-fire risk causing increased disturbance in ecosystem, wildlife, and public community health fare (Cravens et al. 2021;Gupta et al. 2020). Occurrence of drought has been initially attributed to the abnormal meteorological water conditions, however, presently also affected by human water usage managements. Thereby, indirectly affecting the overall agricultural, industrial, and energy sector (Cravens et al. 2021). Further, owing to the fluctuating climatic adversities, drought stress has become a serious threat to plant growth and productivity (Begum et al. 2019;Kaya et al. 2019). Physiological, biochemical, and metabolic alterations, such as overall reduction in photosynthetic efficiency, stomatal conductance, production of essential metabolites involved in the cellular pathways (Farooq et al. 2020;Jha et al. 2020;Kaya et al. 2020) and interference in the symbiotic association of rhizobium with legumes, subsequently affecting root-nodule development and nitrogen fixation are some of the damaging aspects of drought stress in plants, often leading to reduced plant vigour and grain yield (Jha et al. 2020). Moreover, impact of drought stress in the annual loss of agricultural crop yield surpasses the combined effect of combined biotic stress elements, such as pathogen attack (Gupta et al. 2020). Hence, highly climatic resilient crops are the need of the hour for augmenting the agricultural production and profitable agriculture that would enable feeding the ever-increasing global population with healthy and nutritious diet.
Cowpea (Vigna unguiculata L. Walp) is a highly nutritious grain legume crop (diploid; 2n = 2x = 22) and serves as the primary protein source for millions of people in tropical and sub-tropical regions of Africa, Asia, and other developing countries of the world. Due to its ability to associate with nitrogen-fixing bacteria and vesicular-arbuscular mycorrhizal fungi, cowpea can survive in low fertile soils, and often grown with cereals in intercropping system. It is a promising crop for cultivation under adverse climatic fluctuations, being highly resilient to heat and drought stress in comparison to other members of the Leguminosae family (Carvalho et al. 2017). Cowpea germplasm also comprises of sensitive varieties showing differential response under drought stress, and during the vegetative-stage, imposed drought stress affects the overall growth, gas-exchange, stomatal conductance, and the nutrient absorption ability. Therefore, improvement of the cowpea population would be of great importance and further increase its yield in the diverse agro-ecological environments (Nkomo et al. 2021). The mechanism of drought tolerance in cowpea has been a complex trait with the involvement of several genes (Carvalho et al. 2017;Ravelombola et al. 2021). Further, recent advances in high-throughput sequencing have paved the way for transcriptomic study in legumes, like alfalfa Wan et al. 2021), chickpea (Badhan et al. 2018), peanut (Brasileiro et al. 2015), lentil (Hosseini et al. 2021), and soybean (Tamang et al. 2021) to facilitate identification of drought-related gene networks for understanding varied drought response mechanisms in legume crops. Thus, the identification of the differentially expressed drought-responsive regulatory genes and their networks can be achieved by emphasizing on the large-scale sequencing coverage study of cowpea transcriptome.
The microRNAs (miRNAs) comprise of a class of endogenous small non-coding RNAs varying in size from 21 to 24 nucleotide (nt), which are highly conserved across the plant species, and known to regulate plant growth, development, generate signalling cascade for genome-wide response to various biotic and abiotic stresses (Pagano et al. 2021). Further, with the advent of high-throughput sequencing, the existing repository for miRNAs, the potential regulators of gene expression, mainly via post-transcriptional gene silencing, have been expanded across the plant kingdom. The involvement of miRNAs in response to varied environmental cues like salinity, drought, cold, and nutrient deficiency have been reported in several plants, like Arabidopsis (Fasani et al. 2021), tea , potato (Liao et al. 2021), sunflower (Liang et al. 2020), including legumes such as, grass pea (Bhat et al. 2020), pigeon pea (Buch et al. 2020), and peanut (Ren et al. 2020). So far cowpea is concerned, only Barrera-Figueroa et al. (2011) have performed miRNA profiling in response to drought stress in the cultivars contrast for drought tolerance. In addition, Shui et al. (2013) identified miRNAs in cowpea using in-silico approach and performed expression analysis of their targets. However, a comparative study focusing on the tissue-specific expression of miRNAs integrated with transcriptome profiling in response to drought stress is yet to be reported in cowpea. Furthermore, owing to its small genome size (~ 620 Mb) in comparison with the other legumes, cowpea is an ideal model crop system for functional genomic studies. Therefore, our study focused on an integrated approach towards deciphering the drought-induced alterations in the miRNA and mRNA transcriptome profiling of a high-yielding Indian droughttolerant cowpea cultivar, Pusa Komal, and further identified the potential miRNA(s) and corresponding target gene(s) involved in drought tolerance mechanism in the plant.
In this study, we found the differential expression of a total of 198 highly conserved, 21 legume-specific, 14 less-conserved, and 10 novel miRNAs combined in both leaf and root libraries of cowpea. Furthermore, a total of 3391 (up-regulated) and 3799 (down-regulated) differentially expressed genes (DEGs) combining the leaf and root libraries of cowpea were identified in response to drought stress. The most noticeable observation was a very high up-regulation of vun-miR408 in response to the drought stress limited to only the leaves of cowpea, which was accompanied by a significant decrease in the expression of the target gene laccase 12. Furthermore, overexpression of the artificial miR408 in the plant resulted in its increased drought tolerance, suggesting the potential use of such biotechnological intervention in improving drought tolerance of the cowpea cultivars of interest. Taken together, our work provided an insight into cowpea drought tolerance mechanism through omics and translational approach.

Plant material and stress treatment
Healthy cowpea (Vigna unguiculata L. cv. Pusa Komal) seeds were planted in pots filled with soil:sand:compost (2:1:1) mixture and grown in greenhouse under controlled growth conditions (26 ± 2 °C/ 60-70% RH,16 h/8 h light/ dark photoperiod). 1-month-old cowpea plants were divided into two experimental sets and six replicates were kept for each set. One set comprised of control plants that were well-irrigated with leaves showing relative water content (RWC) around 85%. The second set comprised of stress-treated cowpea plants subjected to water-deficit irrigation treatment by withholding watering for 2 weeks and relative water content (RWC) of leaves ranged between 65-66%. Thereafter, the leaf and root tissues from six replicates of control unstressed (US) and drought-stressed (DS) cowpea plants were frozen in liquid nitrogen and stored at -80 °C, until further use.

Differential expression profiling of drought-responsive miRNAs (DE) and target gene prediction
Expression profiling of known and novel miRNAs was performed using mirDeep2 software which generates normalized counts based on information obtained from mapped data. Structural stability and minimum free energy were calculated, which subsequently provided miRDeep scores where a threshold of 0 and above was used as a screening parameter. Differential analysis of miRNA abundance in US and DS leaf and root tissues was conducted using DESeqR package considering p value < 0.05 as a threshold, wherein, DE with log 2 -fold change ≥ 1.0 or ≤ − 1.0 have been considered significantly up-/down-regulated, respectively. The potential targets of DE were obtained using the web-based psRNA Target program (http:// plant grn. noble. org/ psRNA Target/) with default parameters of Schema V2, 2017 release. The selective targets of expectancy 5.0 or less were only considered for annotation and validation.

Comparative expression analysis for selected DE and target genes
For comparative tissue-specific analysis of miRNA expression under drought stress, 40 µg of total RNA isolated from control (unstressed) and drought-stressed leaf (CL, DL) and root (CR, DR) tissues of cowpea was used for northern assay. Probes (DNA oligos) complementary to the selected DE were end-labeled with [γ-32P] dATP (Table S1). Northern blotting and analysis were performed following the method as described by Chandra et al. (2021). First-strand cDNA synthesis was performed using the QuantiTect Reverse Transcription Kit (Qiagen, USA) from the total RNA of US and DS leaf and root tissues of cowpea. Primers listed in Table S1 were designed using Primer Blast software (http:// www. ncbi. nlm. nih. gov/ tools/ primer-blast/). Real-time PCR was performed using QuantiNova SYBR Green RT-PCR Kit (Qiagen, USA) on LightCycler ® 480 Real-time PCR System (Roche Life Sciences, Penzberg, Germany). The qRT-PCR was set using biological duplicates and technical triplicates (n = 6) with β-tubulin gene as an internal control.

3
Validation of miRNA guided cleavage of target genes by 5′ RLM-RACE Using FirstChoice ® RLM-RACE Kit (Thermo Fisher Scientific Inc., MA, USA), four RACE libraries specific for US and DS leaf and root tissues of cowpea were generated. Briefly, 10 µg total RNA from each tissue was ligated to a 5′-RACE adapter (5′-GCU GAU GGC GAU GAA UGA ACA CUG CGU UUG CUG GCU UUG AUG AAA-3′), incubated at 37 °C for 3 h prior to cDNA synthesis, according to the manufacturer's instructions. The nested PCR products generated using RACE Inner and gene-specific (GSP) primers (Table S1) with 1/10 volume of outer race product as template, were subsequently cloned into pGEMT-easy vector, transformed into competent DH5α bacterial cells, and sequenced.

Development of transgenic cowpea lines overexpressing amiR408
The plant binary construct for vun-miR408-3p was obtained by overexpressing amiR408-3p under CaMV 35S promoter in pCAMBIA 2301 (11.6 kb). Briefly, the ath-miR319a/a* sequences in pRS300 vector were replaced with miRNA 5′-ATG CAC TGC CTC TTC CCT GGC-3′ and miRNA* 5′-GCC AGG GAA GAG GCA GTG CAT-3′ sequences using oligos (Table S1) designed by the WMD3 software (http:// wmd3. weige lworld. org/ cgibin/ webapp. cgi/). The aMIR408 precursor fragment was cloned into EcoRI and BamHI sites of pRT101 (3.3 kb). Thereafter, the CaMV 35S promoter-aMIR408-CaMV 35S terminator fragment (~ 1.1 kb) flanked by HindIII sites was further, subcloned into pCAMBIA2301. The resulting pCAMBIA2301::35S-aMIR408 plasmid (~ 12.7 kb) was transformed into Agrobacterium strain EHA105 by freeze-thaw method. The transgenic cowpea lines were generated following Agrobacterium-mediated transformation protocol, previously described by Mishra et al. (2014). The independent T 0 and T 2 generation progenies of OX-amiR408 transgenic lines (1 and 5) were confirmed for the presence of the transgene by southern blotting with nptII and gus-A probe. Briefly, 50 µg of EcoRI digested purified genomic DNA was run on 1% agarose gel and transferred to nylon membrane (Hybond N + , GE Healthcare) by overnight capillary transfer following depurination, denaturation, and neutralization. The 796 bp and 864 bp purified denatured genespecific fragments for nptII and gus-A, respectively, were labeled with [32P] dCTP using the random primer labeling method (Prime-a-Gene Labeling System Kit, Promega, Madison, WI) at 37 °C for 1 h. The membrane was hybridized with nptII specific probe using PerfectHyb Plus hybridization buffer (Sigma-Aldrich, Missouri, US) at 37 °C for 16 h, washed two times at 60 °C for 15 min in non-stringent, 2X SSC + 0.1% SDS and stringent (1X SSC + 0.1% SDS) condition. Thereafter, the membrane was exposed to X-ray film, stored at − 80 °C for 1-2 d prior to being photographed. The blots were stripped and re-blotted using gus-A probe. Similarly, northern blot analysis was performed to detect the expression of amiR408 in southern confirmed transgenic lines. The primers used for transgenic analysis and screening of OX-amiR408 lines are listed in Table S1.

Morphological and physiological changes under water deficit and salt stress
To study phenotypic response to water-deficit stress, cowpea plants were taken in two growth stages. Firstly 14-dayold and secondly 1-month-old untransformed control and T 2 generation progenies of transgenic OX-amiR408 lines 1 and 5 grown in greenhouse under controlled environmental conditions, were subjected to water-deficit irrigation for a period of 15 days and revival was studied after 20 and 15 days, respectively. Further, leaf samples were harvested for total chlorophyll estimation, relative water content, proline content, and detection of hydrogen peroxide (H 2 O 2 ) by 3,3′-diaminobenzidine (DAB) staining, according to Mishra et al. (2014). Further, early-seedlings (3-day-old) of untransformed and T 2 generation of transgenic OX-amiR408 lines 1 and 5 (T 2 .1 and T 2 .5) were grown hydroponically in ½ MS media supplemented with/without 200 mM NaCl for 72 h, and subsequently transferred to fresh media (HiMedia, Mumbai, India), photographed on 7th day of revival. Statistical difference between the variances was estimated by ANOVA (analysis of variance) and statistically significant difference was measured using six replicates (n = 6) by Bonferroni analysis at P ≤ 0.05.

SEM analysis
Leaf samples from untransformed control and T 1 generation progenies of OX-amiR408 lines were used for studying stomatal changes and trichome density, three replicates from each plant were used. Sample preparation was performed by vacuum-infiltrating leaf samples with 2.5% (w/v) glutaraldehyde in 0.1 M phosphate buffer solution (pH 7.4) buffer. After overnight fixation, samples were further dehydrated with a graded ethanol series (50, 70, 80, 90, 100%). Samples were dried in critical-point dryer (Leica Microsystems), mounted on carbon stub and gold-coated before imaging using Scanning Electron Microscope (Zeiss EVO 18 SEM, Oberkochen, Germany) (Bhattacharya et al. 2020).

Transcriptome sequencing, data processing, and expression analysis
Four RNA sequencing libraries (CL, DL, CR, DR) were prepared with Illumina-compatible NEBNext ® Ultra™ II Directional RNA Library Prep Kit (New England BioLabs, MA, USA) according to manufacturer's instructions (Genotypic Technology Pvt. Ltd., Bangalore, India). Similarly, six RNA sequencing libraries prepared from leaf and root of control untransformed wild-type (WT) and T 2 generation of transgenic OX-amiR408 lines 1 and 5 were analysed using Illumina HiSeq platform (Clevergene Biocorp Private Limited, Bangalore, India). Briefly, mRNA isolation was performed from 1 µg total RNA using oligo-dT magnetic beads and further, subjected to fragmentation and priming followed by cDNA synthesis. The cDNA fragments were amplified to generate transcriptome libraries and sequenced on Illumina HiSeq XTen sequencer (Illumina, San Diego, USA) for 150 bp paired-end chemistry following manufacturer's procedure. The sequencing quality was assessed using FastQC v0.11.8 software. Transcriptome analysis was performed by processing the raw data for removal of low-quality data (< q30) and adaptor sequences using TrimGalore (https:// www. bioin forma tics. babra ham. ac. uk/ proje cts/ trim_ galore/). The pre-processed high-quality data were aligned to reference Vigna unguiculata genome using Hisat2 with the default parameters to identify the alignment percentage. HTSeq was used to estimate and calculate transcript abundance and differentially expressed genes (DEG) were calculated using DESeqR package. Transcripts were categorized into up, down and neutrally regulated based on the log 2 -fold change cut-off value of 1. Gene Ontology (GO) annotation was performed by homology search against the Viridiplantae protein sequences from the Uniprot database, with an e-value cut-off of e-5. Pathway analysis was performed using the (KAAS) KEGG Automatic Annotation Server (http:// www. genome. jp/ kegg/). Real-time PCR analysis was performed for the selected DEGs in the untransformed control (wild-type) cowpea plants and T 2 generation progenies of transgenic OX-amiR408 lines 1 and 5 using primers as listed in Table S1. The qRT-PCR was set using biological duplicates and technical triplicates (n = 6) with β-tubulin gene as an internal control.

Transcriptome sequencing and differential gene expression analysis
1-month-old cowpea plants were subjected to natural soil drying process for 2 weeks following which the leaves started wilting and displayed reduced plant growth with a 23% decrease in RWC (Fig. S1). In this study, four libraries (control leaf, CL; control root, CR; drought-stressed leaf, DL; drought-stressed root, DR) were sequenced using the Illumina sequencing platform. A total of 151.77 million paired-end raw data were generated and nearly 96.60% of total reads were retained as high quality (> Q30) and adapter-free processed data. An average of 88.67% of the reads were aligned to the Vigna unguiculata reference genome with 93.35% average alignment for CL and DL libraries and 88.5% average alignment for CR and DR libraries (Fig. S2). Further, an average of 16 590 transcripts were found to be expressed in each of four libraries. Under drought stress, the number of up-regulated differentially expressed genes (DEGs) ranged around 2995/497 in leaf/root tissues, similarly, the number of down-regulated DEGs ranged around 3250/700 in leaf/root tissues as shown in the volcano plots ( Fig. 1a, b, Table S2). There were nearly 101 (up-regulated) and 151 (down-regulated) DEGs common to leaf and root RNA-seq libraries under drought stress, with 20 most significant DEGs in leaf and root shown in the heatmap (Fig. 1c, d). The DEGs were functionally annotated using homology approach (blastX) against Viridiplantae from the Uniport database and gene ontology (GO) terms were assigned. GO enrichment analysis revealed that under biological process category, transcription regulation (GO: 0006355), carbohydrate metabolic process (GO:0005975), signal transduction (GO:0007165), cell wall organization (GO:0071555) and DNA integration (GO:0015074) were the most significantly enriched terms. Under cellular component category, integral component of membrane (GO:0016021), nucleus (GO:0005634), cytoplasm (GO:0005737), and cell (GO:0005623) terms occurred significantly. Under molecular function category, ATP binding (GO:0005524), DNA binding (GO:0003677), metal ion binding (GO:0046872), protein kinase activity (GO:0004672), and sequencespecific DNA binding transcription factor activity (GO:0003700) were the most significantly enriched ones. Further, the pathway analysis for DEGs was done by searching KAAS database. A total of 232 pathways were identified for DEGs and the most enriched occurrence was observed for pathways involving chromosome and associated proteins [ (Fig. 2a, b). Plants display efficient sensory correlations with outside environment to facilitate cellular homeostasis by vital reorganization of solute transport across its membrane using several active and passive transporters under abiotic stress (Conde et al. 2011). Among the droughtresponsive DEGs, several members of solute carrier family proteins, glutathione-S-transferase, aquaporins, major facilitator superfamily (MFS), MATE-type proteins, heat shock proteins (HSP), and PRA-1 family proteins were significantly up-regulated. Similarly, transcription factors previously known as regulators of drought stress signal transduction pathway, constituted 23% in leaf and 26% in root, with a higher occurrence of MYB, NAC, homeobox leucine zipper, heat shock TFs, WRKY, AP2/ERF, EREBP, and NFYA in the leaf and root library in cowpea (Fig. 2c).

Small RNA sequencing and analysis
For identification of drought-responsive miRNAs in cowpea, four small RNA (sRNA) libraries (control leaf, CL; control root, CR; drought-stressed leaf, DL; droughtstressed root, DR) were prepared and sequenced in this study. Overall, 160,981,576 total reads and 20,759,281 unique reads were obtained from all four libraries. After subsequent discarding of nearly 140.22 million reads owing to low quality and size trimming, 5.03 million high quality and non-redundant reads were retained for analysis with 7,925,343 (CL), 7,038,949 (DL), 3,139,238 (CR), and 2,040,451 (DR) reads in each library. Nearly, 90-99% of clean reads aligned to cowpea genome assembly (ASM411807v1). After removal of nearly 1.6 million reads mapping to ribosomal RNAs (rRNAs), small nucleolar RNAs (snoRNAs), small nuclear RNAs (snR-NAs), and transfer RNAs (tRNAs), a total of 12.3 million clustered reads from all four libraries were further used for identification of known and novel miRNAs in cowpea ( Table 1). The percentage of sRNA abundance across the length distribution range of 16-40 nucleotide (nt) in all four libraries, showed 24-nt sRNAs (30.4-56.6%) were more abundant than 23 nt (14.5-47%), 22 nt (18.9-22.6%) and 21 nt sRNAs (13.8-19.03%) (Fig. S3). The presence of high endogenous small interfering RNAs (24-nt siRNAs) in both control and drought-stressed libraries was consistent with previous studies in switch grass (Xie et al. 2014), and foxtail millet ).

Identification of known and novel miRNAs
For identification of known miRNAs in cowpea, clean unique reads from all four libraries were mapped against mature Viridiplantae miRNAs at miRbase 22 database by homology-based approach, with an e-value of e-4 and nongapped alignment. Overall, 250 known miRNAs with 240 (CL), 242 (DL), 222 (CR) and 201 (DR) miRNAs specific to each library were identified (Fig. S4a, Table S3a). In cowpea, 35 highly conserved, 20 legume-specific, and 14 less-conserved miRNA families were detected (Table S3a; Table S5f). Out of 35 highly conserved miRNA families, miR159 family displayed highest members followed by miR171, miR156, and miR396 (Fig. S4b). Again, miR159a and miR1507 were observed to be most abundant in reads. For detection of putative novel miRNAs, the unaligned unique reads were mapped against the cowpea genome with no mismatch and processed using Mireap_0.22b. The prediction was based on alignment, formation of stem loop structure, minimum free energy as calculated with Vienna package for the secondary structure and location of mature miRNA on the precursor arm. The secondary structure was visualized using RNAfold. A total of 837 novel miRNAs were predicted in all four libraries with 198 (CL), 353 (DL), 212 (CR), and 123 (DR) miRNAs in each library (Fig. S4a). The common occurrence of specific known and novel miRNAs in each library was observed (Fig. S4c,  d). In this study, 224 high-confident novel miRNAs with (MFEI ≥ 0.85) and free energies of stem loop structures ranging from − 20 to − 112.2 kcal·mol −1 were obtained 1 3 (Table S3b). Legume-specific miRNAs such as miR4415, miR5037, miR5374, miR5770, miR6300, and miR9726 (Subramanian et al. 2008;Kulcheski et al. 2011;Zhai et al. 2011;Turner et al. 2012;Goettel et al. 2014) reported only in soybean, similarly, mir2199 (Jagadeeswaran et al. 2009) and miR2650 (Lelandais-Brière et al. 2009), so far known in barrel clover, were also detected in all four libraries in cowpea (Table S3a). Owing to the splicing variation, several isoforms of miRNAs are generated and their differential accumulation in distinct tissues of plants have been evident (Srivastava et al. 2015). Similarly, several isoforms (isomiR-NAs), novel variants and novel star miRNAs for known miRNAs were also observed in this study. IsomiRNAs are a result of imprecision or alternate cleavage by Dicer during pre-miRNA processing with additional nucleotide at 5′-or 3′-end of mature miRNA (Kulcheski et al. 2011). IsomiRNAs for miR1507a, miR1511, miR166, miR167, and miR399 were detected in this study. Similarly, novel isoforms for miR172c, miR2111a, miR477c, miR482g, miR4407, miR530b, miR862a, miR1510b, miR2597, miR5225, and miR5261 were identified in the four libraries of cowpea (Table S4a-c).

Target prediction, expression analysis, and validation by 5′ RNA ligase-mediated rapid amplification of cDNA ends (5′ RLM-RACE)
The target prediction for DE in all four libraries was performed using psRNAtarget (Schema V2 release) online tool using reference cowpea cDNA sequences with default parameters (Table S5h). Further, the cleavage of putative targets of selected DE validated in northern analysis were assayed by 5′ RLM-RACE method (Fig. 4b). It was observed that vun-miR167b, vun-miR398c, and vun-miR171m cleaved auxin response factor 8 (ARF8, LOC114188087), copper-chaperone for superoxide dismutase (CCS, LOC114178604) and scarecrow-like protein 6 (SCL 6, Conserved DE members of miR167, miR169, miR171, miR319, miR398, miR408, miR1514, and novel miRNA, vun_miR_27 were detected in control (CL) and drought-stressed (DL) leaf library of cowpea. b Validation of candidate target genes cleaved by selected DE in cowpea by 5′ RLM-RACE. The target genes NAC 29, ARF 8, SCL 6, CCS, LAC 12, and mcfB were obtained as cleaved products of DE vun-miR1514/ vun_novel_796, vun-miR167, vun-miR171, vun-miR408, and vun_novel_792, respectively. The cleavage of the targets was observed at canonical 10/11th position in most cases, but, also at 12/13th position in mcfB. c Real-time PCR analysis of selected target genes confirmed as cleaved products for selected DE in 5′ RLM-RACE approach. The target genes NAC transcription factor 29 (NAC29), auxin response factor 8 (ARF8), copper chaperone for superoxide dismutase (CCS), scarecrow-like protein 6 (SCL6), and mitochondrial substrate carrier family protein B (mcfB) exhibited inverse expressional co-relation with their respective miRNA counterparts, as compared to northern analysis. Data represent relative fold change in expression of candidate genes in the leaf and root tissues of cowpea (mean value of n = 6 ± SE). Statistical analysis was performed using paired t test and significance is indicated by Asterisk (*) at P ≤ 0.05 LOC114166860), respectively, exactly at 10/11th position unlike, vun_novel_792 which cleaved mitochondrial substrate carrier family protein B (mcfB, LOC114194150) at the 12/13th position from 5′-end of miRNA. Again, 5′-cleaved PCR products (2 bands, Fig. S5) for NAC transcription factor 29 (NAC 29, LOC114167348) was obtained at 10/11th position from 5′-end of vun-miR1514a-5p and vun_novel_796. The above-mentioned targets were observed in CL RACE library and their expression levels were high in drought-stressed as compared to control unstressed leaves (Fig. 4c). However, laccase 12 (LAC 12, LOC114177496) was observed to be precisely cleaved at 13 bp upstream of cleavage site of vun-miR408 in DL RACE library (Fig. 4b).

Physiological analysis of transgenic cowpea lines overexpressing amiR408
Transgenic cowpea lines overexpressing amiR408 were generated using cotyledonary nodes as explants with a transformation efficiency of 2% using 300 explants, following Agrobacterium-mediated transformation method as described in Mishra et al. (2014) and shown in Fig. S6. Two independent transgenic OX-amiR408 lines, T 0 .1 and T 0 .5, were confirmed with the presence of single transgene copy, using a 0.796 kb of nptII probe and a 0.864 kb of gus-A probe. Thereafter, the expression of amiR408-3p was observed in T 2 generation progenies of the overexpression lines of T 0 .1 and T 0 .5, using northern blotting (Fig. 5a, b). Moreover, the 3-day-old cowpea seedlings were studied for withstanding salt stress with 200 mM NaCl, as salinity stress is known to be detrimental towards growth at early seedling stage. However, the transgenic cowpea lines survived at 200 mM NaCl/72 h and post-revival to ½ MS media, in contrast to the control cowpea plants which displayed drastic growth defect and no survival (Fig. 5c). Early vegetative stage, 14-day-old cowpea plants were subjected to water deficit in soil for 15 days. It was observed that both control wildtype (WT) and transgenic lines displayed stunted growth under water-deficit condition, however, primary leaves of the wild-type plant, wilted, and yellowed as compared to the transgenic lines (Fig. 5d). Similarly, 1-month-old mature transgenic lines displayed comparatively better survival and faster revival rate with early flowering (Fig. 5e), and also maintained an average of 27.62% and 30.17% higher RWC and chlorophyll content, respectively, as compared to control wild-type plants under drought stress (Fig. 5f). However, under control unstressed condition no statistically significant change in the measurements for RWC and chlorophyll content was observed between the wild-type and transgenic lines (Fig. 5f). The determination of osmo-protectant like proline has been used as a criterion to observe the cellular osmotic adjustment in drought adaptation response in plants (Golldack et al. 2014). Increased proline content under drought stress has been reported in some cowpea cultivars (Hamidou et al. 2007). Similarly, in this study, a 3.74% significant increase in proline content (4.998 µmoles/g FW) in wild-type cowpea plant under drought stress was measured relative to the control unstressed wild-type plants with 1.054 µmoles/g FW proline content. Similarly, an average of 6.3% increase in transgenic lines (6.184 and 6.468 µmoles/g FW: lines 1 and 5, respectively) under drought stress as compared to the control condition of 0.856 and 0.878 µmoles/g FW for transgenic lines 1 and 5, respectively, signified the importance of proline under drought stress. No statistically significant change was observed in the proline content of wild-type and transgenic lines under control unstressed condition (Fig. 5f). Furthermore, DAB assay showed lesser occurrence of H 2 O 2 in the transgenic lines as compared to the control wild-type plants under drought stress (Fig. 5g). The transgenic lines displayed early flowering (Fig. 5e) and a higher seed yield and longer pods (avg. 11-12 seeds/ pod) than control untransformed wild-type (avg. 10 seeds/ pod) under normal condition (Fig. 5h, i) however, no overall significant change was observed in the seed weight (data not shown). The stomatal aperture with an average length of 21.53 and 14.84 µm (Mag. 3.15 K X) with 3 replicates was observed in SEM analysis (Fig. 6) and guard cells displayed higher turgidity in transgenic lines as compared to the wildtype plants (Fig. S7) as observed under light microscope (Leica DM500, Germany).

RNA-seq, and expression analysis of DEGs in transgenic cowpea lines overexpressing amiR408
RNA-seq was performed for identification of differentially expressed genes between untransformed control (wildtype) and T 2 generation progenies of transgenic cowpea lines overexpressing amiR408 and DEGs were identified using DESeq2 package. Further, genes with read count less than 5 were ignored and absolute log 2 fold change ≥ or ≤ 1 and p value 0.05 were considered significant. Out of 28681 tested genes, a total of 5577/2846 and 3120/2042 common DEGs were up-and down-regulated in leaves/roots of control wild-type and OX-amiR408 lines, T 2 .1 and T 2 .5 (Table S6). The expression analysis for targets of amiR408-3p, such as PLANTACYANIN and laccase-like multicopper oxidases laccase 3 (LAC3) and laccase 12 (LAC12), revealed down-regulation in both control untransformed wild-type plants under drought stress and transgenic lines under normal condition (Fig. 7). Among the DEGs, genes known to be positively related to trichome development such as, Transparent Testa Glabra 1 (TTG1), transcription factor Glabra 3-like (GL3) and zinc finger protein Glabrous Inflorescence Stems (GIS3) were found to be up-regulated in the leaves of control (wild-type) cowpea plants under drought and also, in both transgenic lines (Fig. 7, Table S6. 1 3 c, d). In addition, the SEM analysis for control (wild-type) and transgenic OX-amiR408 lines under normal conditions, also revealed an increased trichome occurrence (Fig. 6). The KEGG pathway analysis of the most significant DEGs in transgenic OX-amiR408 lines, revealed transcription factors like, WRKY, ERF, MYB, AP2, bHLH-type, and homeobox leucine zipper and membrane transporters like aquaporins, solute/sugar transporters, and cation/proton exchangers were potentially enriched in transgenic lines (Fig. 8a).
The GO enrichment analysis of transgenic OX-amiR408 lines revealed regulation of transcription, hydrogen peroxide catabolic process, response to oxidative stress, cell wall and T 2 generation progenies of transgenic OX-amiR408 cowpea lines overexpressing amiR408 under drought and salinity stress. a Southern blot of EcoRI digested genomic DNA of control wild-type (C) and transgenic lines 1 and 5, using probes for nptII and gus-A. b Northern blot performed with amiR408 oligo probed with P 32 -γ-ATP to confirm the expression of amiR408 in transgenic lines, and U6 was used as internal control. c At early seedling stage, 3-day-old cowpea seedlings were exposed to salt stress of 200 mM NaCl for 72 h. The revival was scored after 15 days. Scale bar range (top to bottom): 2 cm; 1 cm; 4 cm d Again, 14-day-old cowpea plants were subjected to water deficit in soil for 15 days and revival observed after 20 days. e 1-month-old mature cowpea plants were subjected to water-deficit irrigation for 15 days and revival scored after 15 days of re-watering condition. f Chlorophyll (mg/g. FW), relative water content (%), and proline content (µmoles/g. FW) estimation was performed under unstressed, US (Day 0 stress) and stress treatment, DS (Day 15 stress) in control (wild-type) and transgenic lines. g Detection of hydrogen peroxide (H 2 O 2 ) by 3,3′-diaminobenzidine (DAB) staining in leaves of control untransformed wild-type and transgenic cowpea lines under control unstressed (US) and drought stress (DS). h-i Pods and their respective seeds were observed for morphological difference in control wild-type and transgenic lines under normal condition. Data values indicate mean values with n = 6 ± standard error (SE) and statistical significance was indicated as different letters (a-c) at P ≤ 0.05 using Bonferroni analysis organization, metal ion transport, and metabolic processes were significantly enriched in biological processes category (Fig. 8b). For instance, real-time expression analysis for selected DEGs such as, PER25 (LOC114162544) peroxidase homologue of AT2G41480, involved in lignification process in Arabidopsis (Shigeto et al. 2014(Shigeto et al. , 2015, increased drastically in leaves of control plants under drought and also in transgenic cowpea lines. WRKY75 (LOC114164714, AT5G13080), known to be expressed under Pi deprivation and also, in development of lateral roots in WRKY75 RNAilines in Arabidopsis (Devaiah et al. 2007) was observed to be repressed in transgenic cowpea lines, similar to being down-regulated in salt and osmotic stress in poplar (Zhao et al. 2019). Again, DEGs such as, GDSL esterase/lipase CPRD49 (LOC114183596, AT3G11210), involved in lipid metabolism up-regulated in control plants under drought stress, however, decreased in the transgenic lines. Among NAC transcription factors, NAC62 (LOC114167089), was observed to be up-regulated in control (wild-type) plants under drought stress and also in the transgenic cowpea line 5 plants (Fig. 7, Table S6). Fig. 6 Study of stomata and trichome density in control wild-type (WT) and transgenic cowpea lines overexpressing amiR408. Scanning electron microscope (SEM) analysis revealed higher occurrence of trichomes and reduced stomatal aperture (Mag. 3.15 K X, Scale: 100 µm). The microscopic analysis was performed with three replicates for each plant Fig. 7 Real-time PCR analysis of miR408 target genes in wild-type (WT) and transgenic cowpea lines overexpressing amiR408. The real-time PCR analysis of target genes of miR408, i.e., PLANTA-CYANIN, LAC3, LAC12, and selected DEGs like PER25, WRKY75, CPRD49, and NAC62 was performed in leaf of untransformed control (wild-type) cowpea plants under drought stress (DL) and transgenic lines under normal condition. The relative fold-change in the expres-sion level of the candidate genes were analysed with reference to the untransformed control (wild-type) cowpea plants under unstressed condition (CL). Data representing mean value of n = 6 ± SE. Statistical analysis was performed using paired t test and significance is indicated by different letters (a-d) at P ≤ 0.05 for each gene individually. Segmentation of graph represents each section for individual genes

miRNA-mRNA network in cowpea under drought stress
Drought being a multi-dimensional stress imposes severe phenomenal changes in grain legumes with restricted growth, development, and metabolic functioning. Plant responses to varied environmental cues, instigates complex reprogramming of responsive genes at transcriptional level in the cells. Further, with the advent of highthroughput sequencing techniques, genotype-based/ tissue-specific transcriptome analysis has contributed to the present understanding of molecular patterns, pathways, and processes for adaptive response against resilient environmental factors in plants. Small RNAs are known as potent regulators of plant growth and development under abiotic stress such as salinity, drought, and temperature fluctuations. Variation in tissue/genotype specific miRNA response has been observed under drought stress, owing to molecular plasticity in miRNA abundance and function in several legumes like Phaseolus vulgaris (Arenas-Huertero et al. 2009;Wu et al. 2017), Medicago truncatula (Wang et al. 2011), Lathyrus sativus (Bhat et al. 2020), Glycine max (Kulcheski et al. 2011;Zhou et al. 2020), Cajanus cajan (Buch et al. 2020), and Macrotyloma uniforum (Yasin et al. 2020). In our study, among the significant DE, cowpea isoforms of gma-miR9726, gma-miR169u, decreased 3.66, and 5.87-fold, respectively, and gma-miR5770a increased 11-fold, in leaf library, similar to expression in drought-sensitive and drought-tolerant common bean (Wu et al. 2017). Similarly, pPvu-miR482-5p decreased 8.9-fold in our root library, however, no expression change was observed in common bean under drought stress (Arenas-Huertero et al. 2009). Among the validated stress-responsive miRNAs in our study, miR1514a is an important legume-specific miRNA, particularly induced in the roots of drought-tolerant cultivar of common bean under drought stress and targeted two NAC transcription factors (TF) NAC000 and NAC700, detected by degradome analysis (Sosa-Valencia et al. 2017). Further, it was detected in higher expression in the leaves than roots in northern blots, under control conditions in common bean (Arenas-Huertero et al. 2009), similar to our finding in cowpea. However, miR1514a was detected in CL in northern blot and was found to cleave a NAC TF NAC29, a homologue of NAC000, as detected in 5′ RLM-RACE in CL RACE library in cowpea. This was further confirmed with increased relative expression of NAC29 (7.9-fold) in drought-stressed leaves of cowpea, indicating the potential role of NAC29 for drought tolerance mechanism in cowpea. However, further studies are required to understand the functional regulation of NAC29 in response to drought stress in cowpea. Similarly, miR398 is a universal stress-responsive miRNA observed under oxidative stress, salinity, drought, and abscisic acid stress (Zhu et al. 2011). In cowpea, out of three isoforms of miR398a/b/c, only miR398c was down-regulated in response to drought stress. The major validated targets of miR398 in legumes are cytosolic copper/zinc superoxide dismutase CSD1, chloroplastic CSD2 and copper chaperone for superoxide dismutase (CCS) in Glycine max , mitochondrial cytochrome c oxidase subunit COX5b in Medicago truncatula (Trindade et al. 2010), and protein disulphide isomerase (PDI) in Medicago sativa (Pokoo et al. 2018). Under oxidative stress, SOD plays a vital role against detoxification of reactive oxygen species (ROS) acting as a first line defence system in plants (Mittler 2002;Ahmad et al. 2010Ahmad et al. , 2019Kohli et al. 2019). Moreover, CCS acts as a chaperone which delivers Cu cofactor to activate CSD1 and CSD2 (Chu et al. 2005). In our case, we identified vun-CCS as a cleaved RACE product for miR398c in CL RACE library. miR398c was observed only in CL in northern analysis and vun-CCS up-regulated (8.5-fold) in DL, indicating, cowpea adapts control against oxidative stress as an immediate response mechanism to drought stress. Similarly, gma-CCS was also identified as a confirmed target using 5′ RLM-RACE and transient-GFP dependent gene expression method in Arabidopsis mesophyll protoplast cells for gma-miR398c . miRNA 167 and miR171 were observed to be repressed under drought stress in our study consistent with findings reported in Lathyrus sativus (Bhat et al. 2020) and Medicago truncatula (Wang et al. 2011), respectively. The conserved miR171 family is involved in diverse functions related to plant leaf and root development such as, shoot branching and radial organization of roots, phytochrome signal transduction, lateral organ polarity, meristem formation, vascular development, and abiotic stress response (Zhu et al. 2015;Vakilian et al. 2020). Similar to our findings, miR171targeted SCL transcription factors, the SCL6 homologues involved in nodulation in Lotus japonicus (De Luis et al. 2012) have been confirmed through RACE approach. SCL6 belongs to the GAI, RGA, and SCR (GRAS) family known to be involved in the developmental processes, such as regulation of axillary bud growth. In Arabidopsis thaliana, one of the members of the miR171 family, miR171c negatively regulates the shoot branching by targeting the SCL6 members (SCL6-II, SCL6-III, and SCL6-IV). Again, the reduced shoot branching phenotype is rescued by overexpressing any one of the miR171c-resistant versions of SCL6 members in the transgenic Arabidopsis lines overexpressing MIR171c (Wang et al. 2010). The increased shoot branching is essential in accumulation of nutrients, light interception, crop yield, and overall growth and survival (Barbier et al. 2017). Moreover, in leaves under drought stress indicating its differential role in cowpea shoot development. miR167 is identified as a potential regulator in auxin and ABA signalling under drought stress in maize (Wei et al. 2009), and was observed to target ARF8 in adult leaves, flowers, and inflorescences in Arabidopsis (Wu et al. 2006). Moreover, recently the auxin-responsive factors ARF8 along with ARF6 in regulation of auxin and gibberellin have been shown to be mediating xylem expansion along with suppression of phloem proliferation and differentiation in Arabidopsis (Ben-Targem et al. 2021). Thus, being suggestive of inter-relation of secondary growth development with the involvement of ARF8, and further research can shed light on its role in the drought-regulatory response. Again, miR167-ARF8 module has been suggested to be conserved in plants regulating the adventitious root formation (Cai et al. 2019). We found in our studies, ARF8 was cleaved in control leaves and increased upon drought stress by nearly two-fold in leaves with decrease in miR167 as evident in northern blots, thus suggestive of its importance in drought tolerance regulation in cowpea. Interestingly, a mitochondrial substrate carrier family protein B-like (mcfB), homologue of AT3G55640, was found to be cleaved by 22 nt novel miRNA, vun_novel_792 at 12/13th position in CL RACE library. Further, a 5.8-fold increase in leaf, and a reduction in roots was observed for mcfB under drought stress. MCFs are particularly transporters of wide range of substrates like nucleotides, amino acids, di-and tricarboxylates, cofactors, vitamins, phosphate, and H ions (Nunes-Nesi et al. 2020). According to the expression profile data analysed for Arabidopsis MCFs available at "The Bio-Analytical Resource for Plant Biology database", Nunes-Nesi et al. (2020) showed AT3G55640 to be highly up-regulated under cold, oxidative and salt stress in shoot tissues, thereby suggestive of its importance in abiotic stress.
The multi-stress-response microRNA, miR408 is known to be differentially regulated in several plants like Arabidopsis (Song et al. 2018), wheat (Feng et al. 2013), barley (Kantar et al. 2010), and barrel clover (Trindade et al. 2010). Again, miR408-3p was observed to be up-regulated and more abundant in leaves in our study, in contrast with no change in both drought-sensitive and tolerant cowpea cultivars (Barrera-Figueroa et al. 2011). Further, an abundant novel isoform of miR408-5p (vun_novel_27), also observed in other legumes like chickpea (Srivastava et al. 2015), sacred lotus (Zheng et al. 2013) was detected under drought stress in our study. Till now, LAC12 have been experimentally confirmed as miR408 target in only Arabidopsis ) and grapevine (Leng et al. 2017). However, a cleaved product of LAC 12, 24-nt upstream to 10/11th cleavage site of miR408 was obtained in our study. The overall regulatory function of miR408 was studied by generation of OX-amiR408 transgenic lines in cowpea. The known targets, plantacyanin (LOC114176940), laccase-3 (LOC114193262) and 5′ RLM-RACE identified laccase-12 (LOC114177496) were observed to be down-regulated under drought stress and also in transgenic OX-amiR408 lines. Overexpression of miR408 resulted in enhanced drought tolerance in chickpea (Hajyzadeh et al. 2015), and non-legume crops ryegrass (Hang et al. 2021); moreover, transgenic OX-miR408 Arabidopsis lines exhibited increased drought sensitivity and salinity tolerance (Ma et al. 2015). In cowpea, miR408 enhanced both drought and salinity tolerance in transgenic OX-amiR408 lines. The measurement of drought-related traits, such as relative water content (RWC), proline content, and chlorophyll content was also performed in the transgenic lines. Relative water content is an efficient method to estimate the dehydration tolerance efficiency in plants as it indicates the plant water status maintained under cellular water-deficit condition. RWC has been conveniently used as a screening parameter to effectively differentiate genotypes contrast to different levels of drought tolerance, such as cowpea (Zegaoui et al. 2017), faba bean (Meißner et al. 2021), maize (Badr and Brüggmann 2020), and wheat (Yadav et al. 2019). The transgenic OX-amiR408 lines maintained a higher relative water content as compared to the wild-type plants under drought stress. Similarly higher accumulation of compatible solute, like proline has been observed in drought-tolerant plant varieties such as, canola (Khodabin et al. 2020), safflower (Farooq et al. 2020), and sunflower (Kosar et al. 2021). Further, proline, a stressresponsive versatile plant metabolite, has been demonstrated to be a potential ROS scavenger (Rehman et al. 2021) which might explain for the higher accumulation of proline along with the lesser detection of cellular H 2 O 2 level in the transgenic lines under drought stress. The microRNA miR408 has been reported to up-regulate in response to ABA and salt stress in poplar (Jia et al. 2009), and rice (Mutum et al. 2013). Here, in this study, several ABA synthesis and responsive genes were also altered as evident in the RNAseq of OX-amiR408 lines indicating plausibly, miR408 and ABA are interlinked in cowpea.
Trichomes are specialized epidermal structures known to be involved in drought tolerance in Caragana korshinskii by reducing water loss and thereby, transpiration rate (Ning et al. 2016). Surprisingly, ryegrass transgenic lines overexpressing Os-miR408 (Hang et al. 2021) also displayed increased bristle like trichomes on leaf surface. The protection of plants against UV damage, avoidance of high temperatures, reduction in water evaporation with enhanced water storage performance was suggested towards the enhanced drought tolerance abilities in tomato accessions with many non-glandular trichomes . The positive regulators of trichome formation includes a WD-40 family protein TTG1 and the C2H2 zinc finger family protein GIS3 which controls trichome initiation in Arabidopsis by integrating of hormone signals (gibberellin and cytokinin) Wang et al. 2021). Moreover, increased trichome density observed as an enhanced morphological feature in the leaves of the transgenic OX-amiR408 cowpea lines in the SEM analysis, further taking into consideration the real-time data analysis showing higher expression levels of TTG1 and GIS3, is suggestive for a plausible connection of the maintenance of enhanced cellular water content with higher trichome density and drought tolerance in the transgenic cowpea lines mediated by the overexpression of amiR408.
The DEGs, such as PER25 (also known as Prx25), CPRD49, WRKY75 related to lignification, lipid metabolism, senescence and lateral root development were also studied for differential change in expression under drought stress. Lignin, a phenolic polymer is essential for the formation of secondary cell walls of plant cells and increasing the mechanical strength and hydrophobicity, required in promoting the long-distance water transportation in plants. Lignin polymerization is promoted by the oxidizing of monolignols by enzymes, namely laccases (LACs) and peroxidases (PRXs) (Hoffmann et al. 2020). Further, the reduction of lignin biosynthesis can lead to reduced plant growth, development, and, moreover, reduced efficacy against environmental stress conditions, therefore, a ligninrich cell wall, the first barrier to any external stimuli can be protective against the imposed oxidative stress, an outcome of severe water stress in plant cells (Liu et al. 2018). The peroxidase 25 (PER25) increased significantly in the transgenic lines, thus further extensive work involving the effect on lignification and drought tolerance regulation in cowpea is required to decipher the exact function of PER25. GDSLtype esterase/lipase proteins (GELPs) belonging to the SGNH hydrolase superfamily are known to play crucial roles in plant growth and development, dehydration response in cowpea (Iuchi et al. 1996), and related to amino acid transport and metabolism pathways in Chinese cabbage (Dai et al. 2021). Interestingly, in our study, CPRD49 was observed to be up-regulated in the wild-type cowpea plants under drought stress in contrast to the down-regulation observed in the transgenic lines. This might be attributed to the lower oxidative stress condition maintained in the transgenic lines, however, further work focusing on cowpea CPRDs and their inter-regulatory role with miRNAs and drought response phenomena can shed light on the exact role played by CPRD49 in cowpea. The WRKY gene family includes the plant-specific transcription factors (TF) involved in several abiotic stress responses such as drought, salinity, alkali, temperature fluctuation, and ultraviolet (UV) radiation . Reports on a WRKY TF, WRKY75 indicate its positive regulation of leaf senescence in Arabidopsis thaliana. WRKY75 is favourably induced by age factor, salicylic acid, H 2 O 2 , and several plant hormones. Further, its absence delayed age-dependent leaf senescence and is also known to suppress H 2 O 2 scavenging in Arabidopsis (Guo et al. 2017), which can be further correlated to the decreased expression of WRKY75 in the transgenic cowpea lines in our study. Thus, a correlation between the expressional changes observed in the transgenic lines provide a framework for the association of the differentially expressed genes involved in the metabolic pathways and hormonal interactions. Thereby, suggesting them as the requisite parameters in the maintenance of coordinated drought tolerance network in the candidate legume crop, cowpea.

Conclusions
The study focused on the differential expression analysis of miRNAs and their target genes in the leaf and root tissues of cowpea under drought stress. The miRNAs, vun-miR169b, vun-miR319f, vun-miR408, and vun_novel_27, an isoform of miR408-5p up-regulated, in contrast to vun-miR167, vun-miR171m, vun-miR398c, and vun-miR1514a, observed to down-regulate under drought stress in the leaves of cowpea, as detected by sRNA sequencing and validated by northern analysis. Being a universal multi-stress-responsive miRNA, miR408 was highly abundant under drought stress, particularly in leaf tissue of cowpea. Hence, the study further focused on functional characterization of vun-miR408-3p, by generating transgenic lines by overexpressing vun-amiR408-3p. The transgenic lines were tolerant to drought and salinity stress exhibiting higher relative water content, proline content, lesser H 2 O 2 production, higher trichome density for reduced water loss, and faster recovery rate post-drought stress. The target genes of miR408, mainly, laccases were found to be repressed. Furthermore, RNAseq of the transgenic lines indicated the significant DEGs were involved in the lipid metabolism process. Laccases have been previously reported as glycosylated multicopper oxidoreductases involved in the oxidization of monolignols to facilitate lignin polymerization process. Alterations in lignin composition has been observed in response to varied environmental stimuli in plants; a reduction in lignin biosynthesis was observed in maize under drought stress (Alvarez et al. 2008). Reduction in the level of LAC3 and LAC12 might be contributing to reduced lignin content in cowpea and indirectly serving as an essential trait in cowpea under drought stress. Moreover, engineered Arabidopsis lines with low lignin and xylan have been shown to display improved drought tolerance, accumulated higher ABA content, stomatal closure and reduced water loss under drought stress (Yan et al. 2018). Again, plants with low lignin and xylan content serve as cost-effective feedstocks for biofuel production.
Future studies focusing on functional characterization of the genes involved in the pathways regulating lignin metabolism in cowpea under drought stress are necessary to extend our present findings on drought tolerance mechanisms in the plant. Deciphering the functional role of laccases, particularly LAC12 in cowpea can provide an insight into how maintaining the level of lignification is an essential parameter related to drought tolerance in cowpea.