Transcriptome sequencing and Differential Gene Expression Analysis
One month old cowpea plants were subjected to natural soil drying process for two weeks following which the leaves started wilting and displayed reduced plant growth with a 23 % decrease in RWC (Fig. 1). 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. S1). 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 regulated transcripts (DGEs) ranged around 2995/497 in leaf/root tissues, similarly, the number of down-regulated DGEs ranged around 3250/700 in leaf/root tissues as shown in the volcano plots (Fig. 2, Table S2). There were nearly 101 (up-regulated) and 151 (down-regulated) DGEs common to leaf and root RNA-seq libraries under drought stress. The DGEs 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 sequence-specific DNA binding transcription factor activity (GO:0003700) were the most significantly enriched ones. Further, the pathway analysis for DGEs was done by searching KAAS database. A total of 232 pathways were identified for DGEs and the most enriched occurrence was observed for pathways involving chromosome and associated proteins , transporters , transcription factors , exosome , membrane trafficking  and chaperones and folding catalysts  (Fig. 3a, 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 drought-responsive DGEs, 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 leaf and root library in cowpea (Fig. 3c). Common DGE reported in cowpea in previous reports were also detected in our leaf library. 9-cis-epoxycarotenoid dioxygenase (VunNCED1) and zeaxanthin epoxidase (VunABA1) involved in ABA biosynthesis, cowpea clones responsive to dehydrationCPRD8, CPRD 22 (Iuchi et al. 1996, 2000) were found to significantly up-regulated under drought stress in cowpea (Table S2).
Small RNA sequencing and analysis
For identification of drought-responsive miRNAs in cowpea, four sRNA libraries (Control Leaf, CL, Control Root, CR, Drought Stressed Leaf, DL, Drought Stressed 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 (snRNAs), and transfer RNAs (tRNAs), a total of 12.3 million clustered reads from all four libraries were further used for identification of conserved and novel miRNAs in cowpea (Table 1). The percentage of sRNA abundance across the length distribution range of 16-40 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. S2). 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 (Wang et al. 2016).
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 non-gapped alignment. Overall, 250 known miRNAs with 240 (CL), 242 (DL), 222 (CR) and 201 (DR) miRNAs specific to each library were identified (Fig. S3a, Table S3). Out of 48 conserved miRNA families, miR159 family displayed highest members followed by miR171, miR156, and miR396 (Fig. S3b). 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 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. S3a). The common occurrence of specific known and novel miRNAs in each library was observed (Fig. S3c, 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 (Table S3). 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 S3). 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 (isomiRNAs), 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 S4).
Differentially expressed miRNAs (DE)
The miRNA abundance was represented as absolute read count and normalized counts (Reads per million mapped miRNA reads, RPM = (Absolute count * 1,000,000)/ library size) for known and novel miRNAs in four libraries. The differential miRNAs (DE) in control and drought stress leaf and root tissues of cowpea libraries were generated using DESeq tool considering size factor. Based on normalized read counts, miRNAs with log2fold change (≥ or ≤ 1) and p-value (≤ 0.05), were considered to be differentially expressed (Table S5). Among, 144 DE in cowpea leaves, 81 up-regulated and 63 down-regulated (Fig. 4a). It was observed that miR393a-5p, miR395i, miR396a-3p, miR5770a were highly induced, contrastingly, miR156a, miR169s, and miR2111, miR482a-3p were highly repressed under drought stress (Fig. 4b, Table S5). Similarly, among 89 DE in roots, 42 up-regulated and 47 down-regulated. Further, DE such as miR166a/b/k-3p, miR319b-3p, miR396e, and miR397a were highly induced, however, miR167d, miR171m/q, miR319d/g, miR396a-5p, miR5374-5p and miR6478 were highly repressed under drought stress (Fig. 4b, Table S5). Further, we observed that miR164e-5p, miR166b, miR169b-5p, miR171p/n, miR403 and miR408 were highly up-regulated whereas, miR164g-5p, miR166e-3p, miR319d/g, and miR396a were highly down-regulated among 54 common DE shared in leaf and root sRNA libraries (Fig. 4c, Table S5). 10 novel DE were differentially expressed in all four libraries under drought stress in cowpea. Selected conserved DE validated by northern blotting showed similar expression pattern to sequencing data except miR171m and vun_novel_27, a novel star of miR408 (Fig. 5a, 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 (http://plantgrn.noble.org/psRNATarget/) using reference cowpea cDNA sequences with default parameters (Table S5). Further, the cleavage of putative targets of selected DE validated in northern analysis were assayed by 5’ RLM-RACE method (Fig. 6). 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, LOC114166860) exactly at 10/11th position unlike, vun_novel_792 which cleaved mitochondrial substrate carrier family protein B (mcfB, LOC114194150) at 12/13th position from 5’- end of miRNA. Again, 5’-cleaved PCR products (2 bands, Fig. S4) 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 leaves (Fig. 5b). 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.
Physiological analysis of transgenic cowpea lines overexpressing amiR408
Transgenic cowpea lines overexpressing amiR408 were generated by 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. Two independent OX-amiR408 lines, T0.1 and T0.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 T2 generation overexpression lines of T0.1 and T0.5, using northern blotting (Fig. 6). Early vegetative stage, 14 days old cowpea plants were subjected to water deficit and salinity stress in soil for 2 weeks. It was observed that both control and transgenic lines displayed stunted growth under water deficit condition, however, primary leaves of control plant, wilted and yellowed as compared with transgenic lines (Fig. 7c). Moreover, the 3 days 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 control cowpea plants which displayed drastic growth defect and no survival (Fig. 7a, b). Similarly, one-month old transgenic lines displayed comparatively better survival and faster revival rate with early flowering (Fig. 8a), and also maintained a higher RWC and chlorophyll content under drought stress as compared to control plants under drought stress (Fig. 8b). 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, 21.1 % increase (4.998 µmoles/g FW) in cowpea control and an average of 13.6 % increase in transgenic lines (6.184 and 6.468 µmoles/g FW: line 1 and 5, respectively) signified the importance of proline under drought stress. Furthermore, DAB assay showed lesser occurrence of H2O2 under stress. The transgenic lines displayed early flowering and a higher seed yield (avg. 11-12 seeds/ pod) than control (avg. 10 seeds/ pod) under control condition (Fig. 8c) however, with no overall significant change was observed in 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 and guard cells displayed higher turgidity in transgenic lines as observed under light microscope (Leica DM500, Germany) (Fig. 9).
RNA-Seq, DGE and expression analysis for transgenic cowpea lines overexpressing amiR408
RNA-seq was performed for identification of differentially expressed genes between untransformed control and T2 transgenic cowpea lines overexpressing amiR408 and DGE were identified using DESeq2 package. Further, genes with read count less than 5 were ignored and absolute log2fold change ≥ or ≤ 1 and p-value 0.05 were considered significant. Out of 28 681 tested genes, a total of 5577/2846 and 3120/2042 common DGEs were up- and down-regulated in leaves/roots of control and OX-amiR408 lines, T2.1 and T2.5 (Table S6). The expression analysis for targets of amiR408-3p, such as plantacyanin and laccase-like multicopper oxidasesLAC3 and LAC12, revealed down-regulation in both control untransformed plants under drought stress and transgenic lines (Fig. 10). Among the DGEs, 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 leaves of CK under drought and both transgenic lines (Fig. 10, Table S6. C-D). In addition, the SEM analysis for control and transgenic OX-amiR408 lines under normal conditions, also revealed an increased trichome occurrence (Fig. 9b). The KEGG pathway analysis of the most significant DGE 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. 11a, b). The GO enrichment analysis of transgenic OX-amiR408 lines revealed regulation of transcription, hydrogen peroxide catabolic process, response to oxidative stress, cell wall organization, metal ion transport and metabolic processes were significantly enriched in biological processes category. For instance, real-time expression analysis for selected DGE such as, PER25 (LOC114162544) peroxidase homologue of AT2G41480, involved in lignification process in Arabidopsis (Shigeto et al. 2014, 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 RNAi-lines (Devaiah et al. 2007) in Arabidopsis, and observed to be down-regulated in salt and osmotic stress in poplar (Zhao et al. 2019), was also observed to be repressed in transgenic cowpea lines. Again, DGE such as, GDSL esterase/lipase CPRD49 (LOC114183596, AT3G11210), involved in lipid metabolism up-regulated in control plants under drought stress, however decreased in transgenic lines. NAC62 (LOC114167089), was observed to be up-regulated in control plants under stress and also in transgenic cowpea line 5 plants (Fig. 10, Table S6).