The Cowpea Kinome: Genomic and Transcriptomic Analysis Under Biotic and Abiotic Stresses

The present work represents a pioneering effort, being the first to analyze genomic and transcriptomic data from Vigna unguiculata (cowpea) kinases. We evaluated the cowpea kinome considering its genome-wide distribution and structural characteristics (at the gene and protein levels), sequence evolution, conservation among Viridiplantae species, and gene expression in three cowpea genotypes under different stress situations, including biotic (injury followed by virus inoculation—CABMV or CPSMV) and abiotic (root dehydration). The structural features of cowpea kinases (VuPKs) indicated that 1,293 bona fide VuPKs covered 20 groups and 118 different families. The RLK-Pelle was the largest group, with 908 members. Insights on the mechanisms of VuPK genomic expansion and conservation among Viridiplantae species indicated dispersed and tandem duplications as major forces for VuPKs’ distribution pattern and high orthology indexes and synteny with other legume species, respectively. Ka/Ks ratios showed that almost all (91%) of the tandem duplication events were under purifying selection. Candidate cis-regulatory elements were associated with different transcription factors (TFs) in the promoter regions of the RLK-Pelle group. C2H2 TFs were closely associated with the promoter regions of almost all scrutinized families for the mentioned group. At the transcriptional level, it was suggested that VuPK up-regulation was stress, genotype, or tissue dependent (or a combination of them). The most prominent families in responding (up-regulation) to all the analyzed stresses were RLK-Pelle_DLSV and CAMK_CAMKL-CHK1. Concerning root dehydration, it was suggested that the up-regulated VuPKs are associated with ABA hormone signaling, auxin hormone transport, and potassium ion metabolism. Additionally, up-regulated VuPKs under root dehydration potentially assist in a critical physiological strategy of the studied cowpea genotype in this assay, with activation of defense mechanisms against biotic stress while responding to root dehydration. This study provides the foundation for further studies on the evolution and molecular function of VuPKs.

The present work represents a pioneering effort, being the first to analyze genomic and transcriptomic data from Vigna unguiculata (cowpea) kinases. We evaluated the cowpea kinome considering its genome-wide distribution and structural characteristics (at the gene and protein levels), sequence evolution, conservation among Viridiplantae species, and gene expression in three cowpea genotypes under different stress situations, including biotic (injury followed by virus inoculation-CABMV or CPSMV) and abiotic (root dehydration). The structural features of cowpea kinases (VuPKs) indicated that 1,293 bona fide VuPKs covered 20 groups and 118 different families. The RLK-Pelle was the largest group, with 908 members. Insights on the mechanisms of VuPK genomic expansion and conservation among Viridiplantae species indicated dispersed and tandem duplications as major forces for VuPKs' distribution pattern and high orthology indexes and synteny with other legume species, respectively. K a /K s ratios showed that almost all (91%) of the tandem duplication events were under purifying selection. Candidate cis-regulatory elements were associated with different transcription factors (TFs) in the promoter regions of the RLK-Pelle group. C2H2 TFs were closely associated with the promoter regions of almost all scrutinized families for the mentioned group. At the transcriptional level, it was suggested that VuPK up-regulation was stress, genotype, or tissue dependent (or a combination of them). The most prominent families in responding (up-regulation) to all the analyzed stresses were RLK-Pelle_DLSV and CAMK_CAMKL-CHK1. Concerning root dehydration, it was suggested that the upregulated VuPKs are associated with ABA hormone signaling, auxin hormone transport, and potassium ion metabolism. Additionally, up-regulated VuPKs under root dehydration potentially assist in a critical physiological strategy of the studied cowpea genotype in this assay, with activation of defense mechanisms against biotic stress while responding to root dehydration. This study provides the foundation for further studies on the evolution and molecular function of VuPKs.

INTRODUCTION
Protein kinases (PKs) comprise the largest family of "molecular switches" responsible for orchestrating protein activities. These actors are associated with almost all essential cellular functions, modifying other proteins by chemically adding phosphate groups (i.e., phosphorylation), resulting in their activation (Rauch et al., 2011). Members of the PK superfamily are involved in the signaling pathways of plants in response to both abiotic (Jaggi, 2018) and biotic (Zhu, 2016) stresses, developmental processes (McCready et al., 2020), and others, suggesting their critical participation in plant homeostasis and adaptation responses.
So far, only preliminary studies with less sensitive technologies (SuperSAGE tags) have been developed, addressing the PK gene expression in cowpea (Vigna unguiculata) (Kido et al., 2011). The mentioned species is a major legume commonly grown on poor soils in warm, dry areas of Africa, Asia, South America, and the United States. Its social role is important, being a relevant source of proteins and minerals for millions of people in sub-Saharan Africa and other developing regions (Boukar et al., 2016). Due to its relevance, significant efforts in the omics field have been carried out for cowpea. Recently, its reference genome was published (Lonardi et al., 2019). The mentioned structure has been fundamental for understanding the physiological superiority of this crop under stresses, for mining of genes with biotechnological potential, and for comparative legume genomics studies.
In Brazil, country responsible for 12% of the world's cowpea production (Boukar et al., 2016), the Cowpea Genomics Consortium (CpGC) was established. This biotechnological network sequenced genomes from different genotypes, which are in an active process of assembly. Besides this, the CpGC has assembled transcriptomes (RNA-Seq) for tolerant/resistant cultivars under root dehydration (a critical component of the complex effects of drought stress) and under combined stresses [mechanical injury of leaves followed by viral inoculation of cowpea aphid-borne mosaic virus (CABMV) or cowpea severe mosaic virus (CPSMV)-two important pathogenic viruses of the mentioned crop]. The available transcriptomic data at CpGC and the diversity of treatments and tissues studied, added to the public data of the cowpea reference genome, represent a rich source of biological information.
The present study aimed to explore and characterize the cowpea's kinome. For this purpose, 1,293 loci encoding cowpea kinases (VuPKs) were identified in the reference genome and categorized into groups and families. Additionally, VuPK genes were structurally characterized; their products were analyzed considering subcellular location and protein properties. Further approaches included the study of VuPK expansion mechanisms in cowpea genome; the determination of associated cis-regulatory elements in VuPK promoter regions as well as VuPK gene orthology within Viridiplantae and synteny analysis among legumes of socioeconomic importance also were carried out. Finally, the CpGC transcriptomics data were scrutinized for VuPK identification and gene expression studies. The present work provides significant advances in the studied theme and brings a starting point for further experimental research in cowpea molecular biology.

Kinase Mining and Identification
The VuPK identification and classification were computationally based on two steps: first, using the iTAK (Zheng et al., 2016) software. This tool has the potential to identify up to 150 kinase families by performing two main actions: (a) Kinase identification: if the analyzed sequences have a significant hit (at least 50% of the Pfam domain model) to the protein kinase domains (PF00069, PF07714, or PF00481) in the Pfam database. (b) Kinase classification: after step "a", the identified VuPKs are classified into gene families by comparing their sequences to a set of Hidden Markov Models (e-value cutoff < 1.0e −5 ) developed by Lehti- Shiu and Shiu (2012).
The second strategy was based on a distance analysis using sequences with kinase domain to confirm the iTAK result. For this purpose, a multiple alignment of candidate VuPKs was generated using the ClustalW (Larkin et al., 2007) and a tree, obtained by neighbor-joining (NJ) method (Saitou and Nei, 1987) with a bootstrap of 1,000 replicates (default parameters).
The annotation of a VuPK was considered accurate (bona fide VuPK) when members of the same family, according to the result of the iTAK software, grouped with their peers in the obtained NJ tree. If the iTAK results and the NJ tree were not convergent, the divergent element was indicated as "unknown" (UNK).
Considering the reference genome, the mentioned steps were carried out against the predicted cowpea proteome (Vunguiculata_469_v1.1.protein_primaryTranscriptOnly.fa; Phytozome database V12). Only the longest peptide sequence for each cowpea gene was used. The same steps were carried out to identify bona fide VuPKs in the cowpea RNA-Seq libraries at CpGC. For gene expression studies, all identified VuPK isoforms were scrutinized.
Multiple Collinearity Scan toolkit (MCScanX; Wang et al., 2012) package-downstream analysis mode: "duplicate_gene_classifier"-was applied to classify the origins of the duplicated VuPK genes. The following procedure was used by MCScanX to assign the duplication mechanisms: (1) All genes were initially classified as "singletons" (i.e., no duplicates within the cowpea genome) and assigned gene ranks according to their order of appearance along chromosomes; (2) BLASTp results were evaluated, and the genes with BLASTP hits to other genes were re-labeled as "dispersed duplicates"; (3) in any BLASTp hit, two genes were re-labeled as "proximal duplicates" if they had a difference of gene rank <20; (4) in any BLASTp hit, two genes were re-labeled as "tandem duplicates" if they had a difference of gene rank = 1; (5) MCScanX was executed. The anchor genes in collinear blocks were re-labeled as "WGD/segmental"; (6) so, if a gene appeared in multiple BLASTp hits, it was assigned a unique class according to the order of priority: whole-genome/segmental > tandem > proximal > dispersed.
Synteny analyses between "cowpea versus soybean" and "cowpea versus common bean" were carried out using the same software (downstream analysis mode: "dual_synteny_plotter") and using default parameters.

Ratio of Synonymous and Non-synonymous Substitutions per Site for Tandem Duplicated Genes
The ratio of non-synonymous (K a ) to synonymous (K s ) substitutions (K a /K s ) was used to determine the selection pressure among tandem duplication events. ClustalW 2.0 (Larkin et al., 2007) software first aligned the full-length coding sequences of tandemly duplicated VuPK genes. Then, the K a and K s rates were calculated with the standard genetic code table by the Nei-Gojobori method (Jukes-Cantor model) implemented on MEGA 7 (Kumar et al., 2016).
The K a /K s ratio is an indicator of the selection history on genes or gene regions: (i) if its value is lower than "1, " the duplicated gene pairs may have evolved from purifying selection (also called negative selection; conserves the amino acid sequence), (ii) if K a /K s equals "1, " which means neutral selection (that had no constraint for sequence divergence); and (iii) when K a /K s is greater than "1, " it means positive selection (that led to different peptides).

Orthology of VuPKs
For orthology detection, information available in the Phytozome 2 was retrieved. In the section PhytoMine (an InterMine interface to data from Phytozome), a query was created in the form of a list with the identities of the genes to be analyzed. Ninety-three 2 https://phytozome.jgi.doe.gov assembled and annotated genomes from 82 Viridiplantae species were searched. The result was generated by InParanoid (Sonnhammer and Östlund, 2015). This software uses pairwise similarity scores, calculated using BLASTp between two complete proteomes, for constructing orthology groups using the similarity criteria described by Ostlund et al. (2010).

Candidate Cis-Regulatory Element Enrichment Analysis
Sequences of promoter regions (within 1.0-kb range) of VuPKencoding genes were obtained from the Phytozome v.12 2 database by Application Programming Interface. The motifs (candidate cis-regulatory elements, CCREs) in each promoter region were searched by the software MEME, v5.0.3 3 (Bailey and Elkan, 1994). For each identified motif, the software reports a corresponding e-value. In the present work, an e-value < 10 −2 was adopted as a cutoff point for the characterization of enriched bona fide CCREs. The maximum number of motifs analyzed in the present study for a single VuPK promoter region was 10, and the extension ranged from six to 50 nt.
The TomTom software, v4.11.2 4 (Gupta et al., 2007), was used coupled to the JASPAR (file: JASPAR2018_CORE_plants_nonredundant) database, aiming to assign a transcription factor (TF) to the bona fide CCREs. Such software aligns the enriched bona fide CCREs identified by MEME against a database (JASPAR) of annotated motifs. These alignments were assessed using the following statistical criteria: p-value (cutoff < 10 −2 ) and q-value [false discovery rate (FDR) cutoff < 10 −2 ]. In this work, the identities of TFs associated with the enriched bona fide CCREs were related to the best hit obtained.

Protein Properties and Subcellular Localization Prediction
The molecular weight (MW) and isoelectric point (pI) for all of the VuPK proteins encoded in the cowpea reference genome (Lonardi et al., 2019) were predicted using the JVir Gel online tool (Hiller et al., 2006). The predicted subcellular localizations of the mentioned VuPKs were determined using the CELLO tool (Yu et al., 2006).

Transcriptomics Approach
Biological Material, Experimental Design, and Stress Application Root dehydration assay Seeds of V. unguiculata cv. Pingo de Ouro (considered tolerant to water deficit and drought; Bastos et al., 2011;Rodrigues et al., 2017) were treated with 0.05% (w/v) Thiram (tetramethylthiuram disulfide) and germinated during 2 days at 25 ± 1 • C and 65 ± 5% of temperature and relative humidity, respectively. The seedlings were transferred to a hydroponic system (Rodrigues et al., 2012) (Supplementary Figure 1A) with aerated pH 6.6 balanced nutrient solution (Hoagland and Arnon, 1950). Plantlets were placed in supports so that the roots of the seedlings were wholly immersed in the solution (Supplementary Figure 1A). The plantlets were grown for 3 weeks (developmental stage V3) in a greenhouse under a natural photoperiod of approximately 13/11 h light/dark cycle, temperature of 30 ± 5 • C and 60 ± 10% relative humidity. After this period, root dehydration treatment was initiated by withdrawing the nutrient solution from treated plants (Supplementary Figure 1A). The roots were collected after 25 min (RD25) and 150 min (RD150) after solution removal (Supplementary Figure 1B). The tissue was immediately frozen in liquid nitrogen and stored at −80 • C until RNA extraction. For each treatment, the respective control plants (Cont.25 and Cont.150 ; Supplementary Figure 1B) were maintained in the nutrient solution and subsequently collected. The experimental design was factorial (cultivar vs. extension of root dehydration period) with three biological replicates (RBs) for each control and treatment implemented (Supplementary Figure 1B). Each RB was composed of two individuals.
The seeds of the cowpea 'Pingo de Ouro' genotype (ID: Pingo_de_Ouro_ 1_2) were formally obtained from the 'Cowpea Germplasm Active Bank' of the Agronomic Institute of Pernambuco (AIP; Recife, Pernambuco, Brazil). Prof. Dr. Antônio Félix da Costa (AIP) kindly carried out the species/genotype identification.

Mechanical injury and virus inoculation assays
The mechanical injury and inoculation with CABMV or CPSMV experiments were carried out under controlled conditions in a greenhouse at the Instituto Agronômico de Pernambuco (IPA; Recife, Pernambuco, Brazil) (Supplementary Figure 2A). For the CABMV assay, the resistant genotype IT85F-2687 (Rocha et al., 1996;Oliveira et al., 2012) was used, whereas for CPSMV, the resistant genotype BR-14 Mulato (Cardoso et al., 1990) was analyzed.
The experimental procedure for both assays was conducted separately and in the same way. Both genotypes were sown and cultivated for 3 weeks (V3 developmental stage) under natural photoperiod and temperature ranging from 28 to 32 • C (Supplementary Figure 2A). Young trifoliolate leaves were mechanically injured with carborundum (silicon carbide) to allow virus penetration. After that, the viral inoculum was applied (Supplementary Figure 2A).
Two collection times were performed after the injury/inoculation for each assay (Supplementary Figure 2B): 60 min and 16 h. Each treatment had its respective absolute or mock control (Supplementary Figure 2B). Tissues were immediately frozen in liquid nitrogen and stored at −80 • C until RNA extraction.
The experimental design was factorial (cultivar versus extension of the post-inoculation period) with three RBs for each control and treatment implemented (Supplementary Figure 2B). Each of these RBs was composed of five plants. Each treatment was carried out in an isolated area to avoid the impact of volatile compounds used by the plants for communication.
The differential gene expression of the two tests involving viruses was composed of a combination of two stresses: mechanical injury and viral inoculation. Plant viruses are unable to initiate an infectious process without assistance (from a vector organism or due to certain agricultural practices, for example) since such organisms are incapable of penetrating the plant cell wall. According to Barna and Király (2004), plant viruses do not have specific cell receptors, unlike bacteriophages and viruses that infect animals. Thus, the set 'mechanical injury and viral inoculation' aims to mimic the infectious process in a natural environment.

RNA-Seq Libraries: Synthesis and Sequencing
Total RNA was isolated using the 'SV Total RNA Isolation System' kit (Promega, United States) following the manufacturer's protocol. Agarose gel (1.5%) and Agilent 2100 Bioanalyzer (Agilent Technologies, EUA) were used to evaluate both the concentration and the quality of total extracted RNA. Only samples with RNA integrity number ≥ 8.0 were sequenced. A 'RNAm TruSeq R Stranded LT-Set A' kit (RS-122-2101) (Illumina, San Diego, CA, United States) was employed in messenger RNA purification and cDNA library construction according to the manufacturer's instructions. Paired-end reads, with 100 pb in length, were generated via the Illumina HiSeq R 2500 system using the following kits: 'HiSeq Rapid PE Cluster Kit v2' (PE-402-4002), 'SBS Kit v2 ' (200 Cycle;, and 'TruSeq R Stranded mRNA LT-Set A' (RS-122-2101). All sequencing steps were performed at the Center for Functional Genomics, University of São Paulo (Piracicaba, Brazil).

RNA-Seq Libraries Assembly and Differential Expression Analysis
The 12 RNA-Seq libraries sequenced for the root dehydration assay were assembled together with the 12 from the duet 'mechanical injury and viral inoculation by CABMV, ' in addition to the other 12 from a combination of 'mechanical injury and inoculation by CPSMV.' Such joint assembly allowed to obtain longer and more robust transcripts. Seventy-two replicates (12 biological for root dehydration + 12 biological for CABMV + 12 biological for CPSMV, with two technical replicates each) were submitted to the pipeline.

Gene Ontology Enrichment Analysis
For the proposed analysis, we used the Network-Based Visualization for Omics (NeVOmics) tool (Zúñiga-León et al., 2018) that identifies statistically over-represented biological process Gene Ontology (GO) terms within a given gene/protein set. NeVOmics adopts GeneMerge statistical algorithm to obtain the over-represented functions or categories in the input list. A hypergeometric distribution (p < 0.05) and an FDR correction (FDR < 0.05) were applied to identify the statistically more represented function annotations. For each treatment studied, the whole 'control vs. treatments' transcriptomes were used as background in the NeVOmics enrichment analyses.
NeVOmics uses an input file in plain text containing a list of genes (KEGG gene ID) or proteins (UniProt Entry ID) obtained by the omics approach. Due to the reduced number (48,364, unreviewed-April| 2020) of proteins annotated for V. unguiculata at the UniProt 5 database, the annotation was performed against Arabidopsis thaliana. This plant species has a much more abundant data amount (163877, unreviewed-April| 2020) at UniProt, besides presenting the best information index (16181, manually annotated-April| 2020) among all plants in the database.
To retrieve the A. thaliana UniProt IDs, initially, we obtained proteins from transcriptomes of the CpGC database translated via TransDecoder 6 , a utility included in Trinity software. Then, a BLASTp was performed (cutoff < e −10 ) between the proteomes of A. thaliana (at UniProt) and cowpea (from CpGC). Only the best hit for each alignment was sent for further analysis.
qPCR: Setup, cDNA Synthesis, Efficiency Analysis, and Relative Expression These activities were performed according to the MIQE guidelines (Bustin et al., 2009). To evaluate the transcriptomic data's accuracy, 10 VuPK transcripts (Supplementary Table 1), with up-regulation indicated in RNA-Seq libraries, were randomly selected for qPCR investigation. Three biological and three technical replicates were used per analyzed treatment, aiming to guarantee the statistical reliability of the process. qPCR reactions were performed in 96-well plates at LineGene 9660 (Bioer) using the SYBR Green detection method.
Aliquots of the same total RNA sample sent for sequencing of the RNA-Seq libraries were employed in this step. Possible contamination with genomic DNA (gDNA) was eliminated via digestion with 'RNase-free DNase I.' RNA quantity and quality were evaluated, respectively, using the NanoDrop ND-1000 UV-vis spectrophotometer (Thermo Fisher Scientific) and agarose gel electrophoresis 1% (p/v) stained with Green Blue (LGC, São Paulo, Brazil). For each sample studied, total RNA (1 µg) was reverse-transcribed, generating cDNA, using the "Improm-IITM Reverse Transcriptional System" kit (Promega) with oligo primers (dT) and following the manufacturer's recommendations. qPCR setup, PCR cycling, amplification efficiency assay, primer pair design, and melting curve analysis (Supplementary Appendix 1) were performed according to Amorim et al. (2018).
"Actin" and "ubiquitin-conjugating enzyme E2 variant 1D" (Supplementary Table 1) were selected as reference genes for qPCR root dehydration assay data normalization. These were indicated from tests performed by Amorim et al. (2018) for the same cultivar, tissue, and condition. "F-box protein" and "Polyubiquitin 10" (Supplementary Table 1) were selected as reference genes for qPCR mechanical injury and viral inoculation data normalization. These were indicated from tests performed by CpGC group (data not published) for the same cultivars, tissues, and conditions. The Rest2009 software (standard mode) was used for relative expression analysis of the target transcripts. Such analysis is based on paired comparisons (of target transcript and reference genes under stress conditions and controls) using randomization and bootstrapping-Pair-wise Fixed Reallocation Randomization Test©; Pfaffl et al., 2002). Hypothesis testing (p < 0.05) was used to determine if differences in the expression of target transcripts under control and treated conditions were significant.

Gene Structure and Protein Property Characterization
For the study of the VuPK genes' structural features, 12 parameters (Supplementary Tables 4, 5) were analyzed. This section will cover the number of kinase domains and the number of introns, considered biologically more informative. The average number of each feature for each family as well as the individual data for each VuPK gene are available at Supplementary Tables 4, 5, respectively.
Considering the cowpea kinase families, the average number of kinase domains per protein varied from 1.0 (RLK-Pelle_Singleton, PEK_PEK, BUB, CMGC_CDKL-Crabout 3.0% of the VuPKs families) to 3.8 (RLK-Pelle_RLCK-XI, AGC_RSK-2-about 2.0% of the VuPK families). About 95% of the different VuPK families have, on average, between 1.9 and 3.6 kinase domains (Supplementary Table 4). Concerning the number of introns, 176 (∼13.5%) VuPK genes, associated with 26 different families, did not show this type of sequence. The remaining (1,117 or ∼86.5%) has one or more introns in its genomic structure. The average number of introns per family varied from zero (RLK-Pelle_LRR-VII-3) to 28 (PEK_GCN2) (Supplementary Table 3 The analysis of intra-family characteristics (for families with more than 10 constituent members) showed that the number of kinase domains presents few variations within each family. However, the number of introns is inconstant and can vary by up to an order of magnitude (ex.: RLK-Pelle_DLSV; Supplementary Table 5). A similar fact can be observed at the group level (ex., TKL, RLK-Pelle, CMGC, CAMK, among others; Supplementary  Table 4): almost total conservation of the average number of domains among the families, with subsequent variation in the number of introns.
Subcellular localization determines the environments in which proteins operate. Such parameter influences protein function by controlling access to and the availability of all types of molecular interaction partners. A total of 1,291 bona fide VuPKs were divided into seven cellular compartments (Supplementary Table 6 Table 6). The other compartments showed a heterogeneous constitution of VuPKs. The CELLO tool did not assign subcellular location only to two VuPKs.
The pIs of the bona fide VuPKs varied from 4.10 to 10.00, with MWs ranging from 19,296.34 to 184,869.01 Da (Supplementary Table 6). VuPKs presented highly variable pIs and MWs in intrafamily terms (Supplementary Table 6).

The Landscape of the Gene Duplications of the Cowpea Kinome
Gene duplication has long been regarded as an important evolutionary force that provides abundant raw materials for genetic novelty and can occur by several mechanisms. Our data indicate that no VuPK groups have been expanded from the whole genome duplication mechanism. Only one VuPK, belonging to the CMCG group, was considered singleton. A total of 140 VuPKs showed proximal duplication, with 132 (94.3%) associated with the RLK-Pelle group (Figure 2).
Furthermore, 116 VuPK tandem duplication events were found in the cowpea genome, covering a total of 309 genes, composing 190 duplicated gene pairs (Supplementary Table 7). These gene pairs belonged to five distinct groups (Figure 2): CAMK, CMCG, RLK-Pelle, STE, and TKL. Tandem duplication events with the highest number of duplicated genes occurred within the RLK-Pelle group (Supplementary Table 7).
Finally, the dispersed duplication mechanism was indicated as the main apparatus for the VuPK expansion in the cowpea genome, being the most abundant for the 20 analyzed groups (Figure 2). A total of 843 VuPK genes were duplicated by such a mechanism. Additionally, it is the only procedure responsible for expanding 13 of the 20 VuPK groups (Figure 2).

Diagnosing the Form of Sequence Evolution on VuPK Tandem Duplicated Genes
The forces that drive natural selection can be inferred from the types of nucleotide substitutions in the coding sequence of genes. An informative parameter of the gene evolution under selection has been the ratio of the rate of non-synonymous substitution (K a -that causes an amino acid change) to the rate of synonymous substitution (K s -that does not cause an amino acid change). After pair-to-pair comparison, the K s /K a values varied between 0.09 and 4.67 (Supplementary Table 8) and the K s /K a means of the 190 pairs of duplicated tandem genes was 0.60 (Supplementary Table 8). Ninety-one percent (173) of these gene pairs had K s /K a < 1, suggesting that these genes are under purifying selection. In another context, it is observed that 9% (17) presented K s /K a > 1, implying that positive selection played an important role in the evolution of this reduced amount of VuPK gene pairs.
To further analyze the aforementioned conservation, a more specific study was carried out by contrasting the cowpea genome against the genomes of the two legume species with higher orthology indexes, aiming at mining kinase-anchoring syntenic blocks. The inter-genomic studies between cowpea and soybean indicated 564 blocks of syntenic genes ( Figure 4A). Of these syntenic blocks, 55 had VuPK-encoding loci ( Figure 5A). For the 'cowpea versus common beans' comparison, 192 synthetic blocks were observed (Figure 4B), a lower number than the previous comparison. However, such blocks showed greater FIGURE 3 | Polar bar chart with 20 species with higher orthology indexes to the bona fide VuPKs identified in the cowpea genome. In rectangles below, are the families covered by these 20 species as well as the number of representatives on the graph above. Corresponding colors between items on the polar bar chart and rectangles perform the association species-family. Pv, Phaseolus vulgaris; Gm, Glycine max; Mt, Medicago truncatula; Cr, Capsella rubella; Cc, Citrus clementina; Gr, Gossypium raimondii; Tc, Theobroma cacao; Sl, Solanum lycopersicum; Eg, Eucalyptus grandis; Me, Manihot esculenta; Ma, Musa acuminata; Kl, Kalanchoe laxiflora; Mg, Mimulus guttatus; Sp, Salix purpurea; Cs, Citrus sinensis; Pt, Populus trichocarpa; Rc, Ricinus communis; Kf, Kalanchoe fedtschenkoi; Tp, Trifolium pratense; Vc, Volvox carteri. length ( Figure 5B). Of these, 17 anchored VuPK-encoding loci ( Figure 5B).

Enriched CCREs Mining in VuPK Gene Promoters: A Focus on RLK-Pelle Group Families
The following factors justify the focus on the RLK-Pelle group:   Qualitatively, the families scrutinized showed some specificities of enriched bona fide CCREs and associated TFs (Figure 6 and Supplementary Table 10). Otherwise, it was observed that enriched bona fide CCREs associated with C2H2 transcription factor class (of which the Dof-type family is part) were present in all families studied, except in RLK-Pelle_CrRLK1L-1 (Figure 6 and Supplementary Table 10).

Transcriptomics of VuPKs Under Stressful Conditions: General Aspects
The metrics of all RNA-Seq libraries are available in the Supplementary Appendix 2.
CpGC RNA-Seq libraries were used to analyze the VuPK expression under unfavorable conditions. Specifically, three distinct conditions were studied: (1) Mechanical injury followed by CABMV virus inoculation.
The numbers of differentially expressed VuPKs for each assay are shown in Figure 7, indicating a substantial discrepancy in the up-regulated VuPK amounts under different conditions/treatments. Considering the combined stresses 'mechanical injury + virus inoculation' in the leaves, VuPKs tend to be up-regulated in the initial times (Figure 7). In turn, root dehydration stress recruited a higher VuPK transcriptional response in cowpea, presenting an increased tendency to an up-regulation at the later treatment time after stress application (Figure 7).

VuPKs Transcriptomics on the "Mechanical Injury and Virus Inoculation" Assays: Crosstalk and Specificities
There were no enriched GO terms for the up-regulated VuPKs set in each 'mechanical injury and virus inoculation' assay. Thus, other strategies were used to characterize and compare these VuPKs. The up-regulated kinome analysis in 'mechanical injury + CABMV inoculation' and 'mechanical injury + CPSMV inoculation' showed us that, despite the presence of some similar families, their up-regulated kinome was peculiar to each assay (Figure 8 and Supplementary Table 11).
The RLK-Pelle_DLSV family (with the largest number of upregulated kinases for 'mechanical injury + CABMV inoculation' and the second for 'mechanical injury + CPSMV inoculation'; Figure 8) represents well the peculiarity mentioned previously. Despite appearing in the transcriptome of up-regulated VuPKs in both conditions mentioned, most up-regulated kinases were restricted to the 'mechanical injury and CABMV inoculation' assay. This family presented crosstalk restricted to two upregulated transcripts (Supplementary Table 11). Another example is the CAMK_CDPK family. Although it appears as up-regulated transcripts in both comparisons, they regard different isoforms (Supplementary Table 11), that is, it did not display common members up-regulated. Overall, only 12 upregulated transcripts, associated with eight VuPK families (i.e., RLK-Pelle_DLSV, RLK-Pelle_RLCK-VI, RLK-Pelle_LRR-Xb-2, STE_STE11, CAMK_CAMKL-CHK1, CMGC_CDK-CRK7-CDK9, RLK-Pelle_LRR-III, and RLK-Pelle_LRR-XII-1) exhibited FIGURE 7 | Amounts of bona fide VuPKs differentially expressed in the studied treatments. VuPKs, Vigna unguiculata protein kinases; MI, mechanical injury; CABMV, cowpea aphid-borne mosaic virus; CPSMV, cowpea severe mosaic virus; RD, root dehydration; min, minutes. Table 11). This number is relatively low considering the set of up-regulated transcripts. Thus, it is suggested that the expression of VuPKs in cowpea was virus dependent and/or genotype dependent since these elements vary between treatments.

VuPK Response to Root Dehydration Stress
Root dehydration stress showed a higher number of differentially expressed VuPKs (Figure 7) compared to the other assays. Thus, GO enrichment analyses were possible, with consequent identification of enriched terms. In the present work, we chose to study enriched GO terms related to the 'biological process' category, aiming to understand the cellular processes activated by the up-regulated VuPKs.
More than three dozen enriched terms were found for T25 (Supplementary Table 12) and for T150 (Supplementary Table 13). A large portion of these covers generic processes traditionally associated with protein kinases, such as 'protein phosphorylation, ' 'signal transduction, ' 'stress-activated protein kinase signaling cascade, ' etc. (Supplementary Tables 12, 13). Such obviousness was eliminated, and an interaction network was adapted, presenting the biologically informative enriched GO terms and the relationship between their constituent elements. For T25, such a network is shown in Figure 9 and in Figure 10 for T150.
The study of biological processes potentially activated by up-regulated VuPKs included important actions implemented during the adaptation process to root dehydration in the tolerant cowpea genotype. In both treatments (T25 and T150), enriched terms associated with the ABA hormone were observed, such as 'positive regulation of abscisic acid-activated signaling pathway' (PA) (Figures 9, 10), indicating a role for the up-regulated VuPKs in this process. Some component elements of the 'PA' term have been associated with heavily modulated transcripts (red dots; Figures 9, 10) as Q38997 (T25; Figure 9), an SNF1related protein kinase, and Q8RWH3 (T150; Figure 10) and YAK1 (a dual-specificity tyrosine-regulated protein kinase). This high modulation was observed for several elements that make up other enriched biological processes (red dots; Figures 9, 10), which have an important impact on plant physiology (see "VuPK Transcriptomics Under Stressful Conditions" in the "Discussion" section).
There were also potentially multi-stress VuPKs. These VuPKs have been characterized as responsive to other stress types, such as those associated with enriched GO terms 'response to cold' (RC), 'response to osmotic stress' (RO), and related to the salt stress response [as 'positive regulation of response to salt stress' (PS) and 'response to salt stress'(RS)] (Figures 9, 10). In both treatments, VuPKs were still associated with terms related to the response to biotic stresses (central nodes with the black border highlighted in Figures 9, 10), including 'defense response' (DR), 'defense response to bacterium' (DB), and 'defense response to oomycetes' (DO).
Specifically for the T25 treatment, VuPKs potentially activated processes related to auxin transport (AT) as well as those related to potassium ion homeostasis (PH) (Figure 9). For T150 treatment, up-regulated VuPKs were involved in signaling processes associated with the Brassinosteroid hormone due to the presence of enriched GO terms such as 'brassinosteroid mediated signaling pathway' (BS) and 'detection of brassinosteroid stimulus' (DBS) (Figure 10).

RNA-Seq Data Validation by qPCR
To analyze the reliability and robustness of in silico gene expression (RNA-Seq libraries), a VuPK sample was also scrutinized by qPCR. Ten up-regulated VuPKs (of which three exhibited crosstalk response among the three performed assays) were randomly chosen to analyze the RNA-Seq data quality. The transcript sample covered eight different kinase families (Supplementary Table 1).
The primer pairs (for target transcripts and reference genes) had an amplification efficiency ranging from 90.4 to 109.8% (Supplementary Table 1). The specificity was confirmed by the FIGURE 9 | Biological process enriched for T25 up-regulated VuPKs showing the interaction network (colored lines) of the associated elements (proteins). The colored dots (green, black, and red) in the protein Uniprot ID correspond to the fold-change log (log 2 FC) found in the present study and are according to the presented scale. The colors of the lines are associated with the colors of the central nodes that contain the abbreviation of the enriched biological process. Different colored lines connected to a given protein (colored dots-red, black, or green) indicate that this protein participates in different enriched processes. Central nodes with the border highlighted in black represent GO terms associated with the biotic stress response. T25, 'RD25 vs. Cont.25'; GO terms enriched (AT, auxin transport; CC, cellular response to chitin; CA, cold acclimation; DR, defense response; DB, defense response to bacterium; DO, defense response to oomycetes; PA, positive regulation of abscisic acid-activated signaling pathway; PS, positive regulation of response to salt stress; PH, potassium ion homeostasis; RC, response to cold; RM, response to mannitol; RO, response to osmotic stress; RW, response to wounding; RD, root dehydration; Cont., control). Adapted from NeVOmics tool output.

Genomic Aspects of the Cowpea Kinome: A General View
The recent publication of the cowpea genome (Lonardi et al., 2019) opened up new possibilities and perspectives for the analysis of cowpea kinases, making the present work a pioneer for the mentioned crop. The presence of 1,293 bona fide VuPKs was found in the cowpea genome, corresponding to about 4% (1,293 out of 31,948) of the cowpea protein-coding loci. Compared to two recently studied kinomes, this value was similar to that found in grapevine (3.7%; Vitaceae; Zhu et al., 2018a) and superior to that reported in pineapple (2.8%; Bromeliaceae; Zhu et al., 2018b). This number is similar to that found in A. thaliana (3.6%; Zulawski et al., 2014). For the Leguminosae family, to which the cowpea belongs, data is available only for soybean, with 4.67% of its protein-coding loci associated with kinases (Liu et al., 2015), a higher proportion than that found in the present study.
FIGURE 10 | Biological process enriched for T150 up-regulated VuPKs showing the interaction network (colored lines) of the associated elements (proteins). The colored dots (green, black, and red) in the protein Uniprot ID correspond to the fold-change log (log 2 FC) found in the present study and are according to the presented scale. The colors of the lines are associated with the colors of the central nodes that contain the abbreviation of the enriched biological process. Different colored lines connected to a given protein (colored dots-red, black, or green) indicate that this protein participates in different enriched processes. Central nodes with the border highlighted in black represent GO terms associated with the biotic stress response. T150, 'RD150 vs. Cont.150'; GO terms enriched (AS, abscisic acid-activated signaling pathway; BS, brassinosteroid mediated signaling pathway; CA, cold acclimation; DR, defense response; DB, defense response to bacterium; DO, defense response to oomycetes; DBS, detection of brassinosteroid stimulus; IR, innate immune response; PA, positive regulation of abscisic acid-activated signaling pathway; RDB, regulation of defense response to bacterium; RC, response to cold; RO, response to osmotic stress; RS, response to salt stress; RD, root dehydration; Cont., control). Adapted from NeVOmics tool output.
All the 20 groups of kinases detectable by iTAK tool were found in cowpea. Of the 150 families, 118 were recognized. RLK-Pelle (about 70% of the analyzed kinome) was the group with the greatest abundance and diversity of families ( Figure 1A). This occurred in several other species, such as soybean (Leguminosae family; Liu et al., 2015), grapevine (Vitaceae; Zhu et al., 2018a), pineapple (Bromeliaceae; Zhu et al., 2018b), and Eucalyptus grandis (Myrtaceae; Lehti-Shiu and Shiu, 2012) among others.
The scientific literature indicates that the high number of RLK-Pelle kinase coding loci in genomes is characteristic of land plants (Lehti-Shiu et al., 2009. The reduced RLK-Pelle sizes in animals and green algae are likely evidence that the group was very small before the chlorophyte lineage diverged from land plants and related charophyte algae (Hedges, 2002). The great group expansion has been postulated to be crucial for plantspecific adaptations because its genes play roles ranging from growth regulation to defense response [for a review, see Boller (2012) and Stahl and Simon (2012)].
Considering the families, RLK-Pelle_DLSV (174) was the most abundant ( Figure 1A). This characteristic (high RLK-Pelle_DLSV abundance) remains at the taxonomic genus or even intra-family levels when analyzing some plant groups. In the first case, it was observed that the RLK-Pelle_DLSV family also contained the largest members among the cotton species (Gossypium raimondii, Gossypium arboretum, Gossypium hirsutum, and Gossypium barbadense) (Yan et al., 2018). In the second case, data indicate that, for the Leguminosae family [e.g., soybean and M. truncatula (Lehti- Shiu and Shiu, 2012;Liu et al., 2015), besides cowpea herein analyzed], RLK-Pelle_DLSV is also the most abundant family. According to Liu et al. (2015), A. thaliana, rice (Oryza sativa), and maize (Z. mays) follow the same abundance tendency. The exceptions regard the grapevine (Zhu et al., 2018a) and pineapple (Zhu et al., 2018b) kinomes.
Considering the spatial distribution, it was observed that, despite the high quantity and diversity of VuPKs on chromosome '3' (Vigun03g; Figure 1B), a random pattern was found, being distributed throughout all chromosomes (and even scaffolds). A similar fact is observed concerning soybean (Liu et al., 2015), grapevine (Zhu et al., 2018a), and pineapple (Zhu et al., 2018b).

VuPKs: From Gene Structure to Protein Properties
We observed that the proportion of intronless VuPK genes (13.5%) was much less than that of total VuPK genes with introns (86.5%), suggesting that intron gain has significantly contributed to the structural divergence of the VuPKs. A similar fact was found in soybean (Liu et al., 2015), a plant from the same cowpea family. However, about ∼ 67% (127/190 comparisons; Supplementary Table 8) of tandem duplicated genes have identical intron numbers. This result indicates that, in cowpea genome, paralogous arising from duplication by other mechanisms may be more subjected to an expansion of these quantitative structures. A possible explanation for this fact would be that tandem duplicated paralogous are located in an identical genomic context. This would be a restrictive factor, maintaining the number of introns. In accordance with the statement above, the high rate of VuPK duplication of the 'dispersed' type (the main apparatus for the VuPK expansion in the cowpea genome; Figure 2) would not suffer the restrictions of a similar genomic context, potentiating the high heterogeneity of the number of introns between VuPK members of the same family. However, such an observation needs further investigation.
Another parameter of biological importance is associated with the number of kinase domains. Usually, domain sites are responsible for a particular function or interaction, contributing to the overall protein role. About 95% of VuPK families have more than one domain in their protein structure. Additionally, at the intra-familial level, domain number and type show little variation. Members in each family generally shared similar conserved domain arrangements, suggesting a common evolutionary history, as reported for pineapple by Zhu et al. (2018b). Ding et al. (2009) observed that kinases could form homo-or heterodimers. Thus, the multiple kinase domain structure of VuPKs may be functionally equivalent to dimerized VuPKs, as suggested by Liu et al. (2015). These authors also suggested that the cooperative activity of multiple kinase domains may be required for specific substrates.
In respect of subcellular location, it was observed that almost all of the VuPKs allocated in the plasma membrane were from the RLK-Pelle group (Supplementary Table 6). These proteins play important roles in perceiving the extracellular ligands and activating the downstream pathways [for a review, see Osakabe et al. (2013)]. They can act as transmembrane proteins, with extracellular receptor domains and intracellular kinase domains (Minkoff et al., 2017). The remaining groups, in addition to some members of the RLK-Pelle group, were distributed throughout the vacuole, chloroplast, mitochondria, cytoplasm, extracellular space, and nucleus. This demonstrates the dynamics and plasticity of these proteins in cowpea, suggesting their action at multiple intracellular sites.
Finally, concerning the VuPK biochemical properties, pI and WM, it was observed that they were variants, even within the family (Supplementary Table 6). This result contrasts with those available for this type of analysis. For instance, in grapevine, Zhu et al. (2018a) observed that such parameters within the same family were generally similar.

The VuPK Gene Duplication Mechanisms in Cowpea Genome
The present study demonstrated that three main mechanisms, with different levels of performance, were responsible for the expansion of this gene superfamily (Figure 2): tandem, proximal, and dispersed gene duplications. Such mechanisms occurred mainly in the RLK-Pelle group, reinforcing the trend found in land plants. As mentioned earlier, this group expanded continuously throughout land plant evolution. Along the evolutionary scale, it is observed that the early diverging bryophyte species, Physcomitrella patens (moss), has 329 RLK-Pelle (Lehti- Shiu et al., 2009). In contrast, there are 610 RLK/Pelle genes in A. thaliana (Shiu and Bleecker, 2001), 872 in grapevine (Zhu et al., 2018a), and 1,418 in soybean (Liu et al., 2015).
The most abundant family in the RLK-Pelle group was the RLK-Pelle_DLSV, which comprises about 20% of that group's members. There were also families with very small numbers, such as RLK-Pelle_C-LEC and RLK-Pelle_Singleton (one locus each; Supplementary Table 4). This indicates that, in cowpea, selective forces act differently on different kinase families. In the legume soybean, the size of the RLK-Pelle subfamilies was also highly variable (Liu et al., 2015).
The dispersed duplication mechanism was mainly responsible for the expansion of the cowpea kinome. In all groups analyzed, the mentioned mechanism was predominant, being the only one responsible for the expansion of 13 specific groups (Figure 2). According to Wang et al. (2012), copies duplicated by dispersion mechanism might arise from transposition, such as 'replicative transposition, ' 'non-replicative transposition, ' or 'conservative transposition.' The dispersed duplicated genes are prevalent in different plant genomes . The third most active duplication mechanism was the 'proximal' type (Figure 2). This was closely associated with the duplication of the RLK-Pelle group. About 94% of such gene duplications occurred in that group. A gene duplicated by this mechanism might arise from small-scale transposition or tandem duplication and insertion of some other genes (Wang et al., 2012).
Both aforementioned duplication mechanisms-responsible for duplicating about 76% of the cowpea kinome-could be influenced by transposable elements. It is known that all families of transposons are capable of duplicating genes or gene fragments, and such events have been detected in a broad spectrum of organisms (Cerbin and Jiang, 2018). About 39% of the cowpea genome comprises transposable elements, which are largely responsible for the size differences of the genomes of species of the genus Vigna (Lonardi et al., 2019).
There was the tandem duplication mechanism also, which occurred by expanding VuPKs from five different groups (CAMK, CMCG, RLK-Pelle, STE, and TKL). This mechanism represented the second largest expansion force of the cowpea kinome, covering 309 VuPKs. However, the largest tandem duplication events involved members of the RLK-Pelle group, both in tandem gene extension and in quantity (it covered 180 of the 190 duplicated gene pairs) (Supplementary Table 8). In cowpea, tandem duplicated genes accounted for ∼24% of the kinome. This percentage is one of the largest reported so far, being higher than that in Arabidopsis (9.5%; Zulawski et al., 2014), soybean (10.6%;Liu et al., 2015), pineapple (12,5%; Zhu et al., 2018b), and maize (17.2%; Wei et al., 2014) but lower than that of grapevine (∼34%; Zhu et al., 2018a).
Finally, the K a /K s substitution rate ratio of the tandem duplicated genes (190 gene pairs) was analyzed. This ratio is an indicator of the selection history of genes or gene regions. The vast majority (91%) of the analyzed gene pairs have K a /K s < 1, that is, they are under purifying selection, which preserves the amino acid sequences. This finding also suggests that duplicated pairs may have some functional redundancy. This mechanism potentially compensates for stochastic loss-offunction in specific gene family members or paralogous in the course of evolution.

VuPK Comparative Genomics: From Punctual Orthology to Syntenic Blocks
A total of 64 out of 93 genomes analyzed presented orthology indexes ranging from 30% (M. polymorpha; Marchantiaceae) to 88% (P. vulgaris; Leguminosae) of the 1,293 bona fide VuPKs (Figure 3). This demonstrates a conservation gradient of VuPKs, opening new perspectives of studies concerning their function and influence on the species in which they were conserved.
Leguminosae members had the three highest orthology indexes (Figure 3), a fact already expected since orthology indexes reflect the phylogenetic distance between species (Kaló et al., 2004). The greatest orthology index occurred with common bean, which is considered a cowpea sister group, and both present some divergence when compared to soybean (Lei et al., 2016).
The analysis of orthologous loci among 'cowpea versus soybean' and 'cowpea versus common bean' showed that orthologous VuPKs participated in a greater number of syntenic blocks with soybean. This was because 'cowpea versus soybean' syntenic blocks were more fragmented, smaller, and more numerous. The comparison 'cowpea versus common bean' revealed longer, more complex, and a lower number of syntenic blocks. According to Ganko et al. (2007), 'dispersed' duplicates, like most of the VuPKs, are paralogous that are neither near each other on chromosomes nor show conserved synteny. However, the predominantly dispersed nature of the VuPK duplication in cowpea genome did not preclude their anchoring in syntenic blocks compared with the other analyzed legumes.

CCRE Enrichment and Associated Transcription Factors on RLK-Pelle Group Gene Promoters
A total of 62 bona fide CCREs associated with 49 different TFs were identified (Figure 6 and Supplementary Table 10). The types and quantity of bona fide CCREs varied in the different families analyzed, indicating that they participate in different processes coordinated by different transcription factors. This is in line with the heterogeneous range of functions reported for kinases: from abiotic stress response (Ye et al., 2017) and developmental processes (Stahl and Simon, 2012) to plant defense (Lehti-Shiu et al., 2009). Interestingly, the class of C2H2 transcription factors was closely associated with families in the RLK-Pelle group, being enriched in eight of the nine families in the analyzed group. These TFs function as key transcriptional regulators in plant responses to a broad spectrum of stress conditions, including extreme temperatures, salinity, drought, oxidative stress, excessive light, and silique shattering among others [for a review, see Wang et al. (2019)]. Otherwise, some members of the C2H2 class, such as the Dof-type (also enriched for some of the families analyzed), are involved in pathogen defense, such as viruses (Kang et al., 2016). In addition to suggesting VuPK participation in the referred processes, such data assist in characterizing the transcriptional dynamic of these proteins under unfavorable conditions.

VuPK Transcriptomics Under Stressful Conditions
The qPCR assay corroborated the in silico RNA-Seq gene expression data for the scrutinized targets. This indicates good practices during the sequencing and synthesis of the analyzed libraries. It also reinforces the rigor of transcriptome assembly and statistical methods of the differential expression analysis. All these elements together ensure the robustness of the information available in the present work.
We propose that the VuPK transcriptional orchestration occurs in a stress-dependent and/or genotype-dependent and/or tissue-dependent manner since the crosstalk was relatively reduced, and the mentioned elements (stress, genotype, and tissue) varied among the three assays (see "Materials and Methods" section). To date, only tissue-dependent and stressdependent kinase expression has been reported for angiosperms. Concerning the first one, Zhu et al. (2018a) reported that the tissue-specific RNA-Seq analysis of pineapple kinase genes showed various expression patterns. Regarding the second case, Zhu et al. (2018b) reported for grapevine that the expression of family members as CMGC_CDK-Pl, RLK-Pelle_LRR-XIIIb, and RLK-Pelle_RLCK-IXa were down-regulated in response to salt, PEG, and drought treatments but that these were up-regulated in response to heat stress; Yan et al. (2018) also reported stressspecific up-regulation for G. hirsutum.
In terms of VuPK families, the expression of RLK-Pelle_DLSV and CAMK_CAMKL-CHK1 stood out. These VuPK families are among the three most up-regulated in all the transcriptomics assays performed in the present work. There are reports of upregulation of CAMK_CAMKL-CHK1 in response to stresses that cause dehydration, such as salt, PEG, and drought (Zhu et al., 2018b). However, there was still no information in the literature about the participation of isoforms of these families in response to the combination 'mechanical injury + virus inoculation'. For RLK-Pelle_DLSV, there were still no reports of its participation in plant stress response (abiotic or biotic). Such results raise certain members of these two families to a higher importance level for functional characterization tests. Additionally, they can become biotechnological targets of great interest in cowpea and related species.
Specifically, in relation to root dehydration stress, 100 families of VuPK had at least one up-regulated isoform throughout the assay. GO terms associated with the positive regulation of the ABA hormone ('PA' and ' AS'), a key hormone in response to drought in plants, were enriched in T25 (Figure 9) and/or T150 (Figure 10). According to Fernando and Schroeder (2016), the role of ABA in drought and salt stress is twofold: water balance and cellular dehydration tolerance. Proteins such as Q38997 (Figure 9), an SNF1-related protein kinase, were among the most modulated for the 'PA' term in T25. Umezawa et al. (2004) indicated that a SNF1-related protein kinase is a positive regulator of drought tolerance in Arabidopsis roots. Concerning T150, the protein Q8RWH3 (Figure 10), a YAK1 (a dual-specificity tyrosine-regulated protein kinases), was one of the most modulated for the term "PA." Results of the work by Kim et al. (2016) indicated that AtYak1 plays a role as a positive regulator in ABA-mediated drought response in Arabidopsis.
It was also observed that GO terms associated with the biotic stress response were enriched in both treatments (central nodes with the border highlighted in black; Figures 9, 10). This may be related to a natural strategy of the Pingo de Ouro cultivar. Prasch and Sonnewald (2013) suggested that abiotic stress factors significantly altered turnip mosaic virus-specific signaling networks, which led to the deactivation of defense responses and a higher susceptibility of Arabidopsis plants. Sinha et al. (2019), in turn, also observed an increased incidence of fungal diseases in chickpea plants under severe drought stress compared to those in well-irrigated field conditions. Thus, the studied cultivar may use defense mechanisms during an unfavorable abiotic condition, aiming to keep on constant alert against different stressors, with VuPKs being important players in the perception and signaling process.
Terms such as 'response to cold' , 'response to osmotic stress' , and those related to salt stress response (RS and PS) (Figures 9, 10) were also enriched in T25 and/or T150. This enrichment was possibly because the mentioned conditions (cold and osmotic and salt stress) strongly alter the water availability in the environment. Chinnusamy et al. (2007) observed that cold temperature affects membrane fluidity and disrupts nucleic acid and protein structures besides hindering water and nutrient uptake. Torres et al. (2007) suggest that osmotic and salt stress should be treated as a group of factors imposing alterations in the plant water status. Thus, the up-regulated VuPKs associated with these terms could act in processes related to signaling and response to water availability. An example is Q39193 (a SnRK2), one of the most modulated for the term 'response to salt stress' (RS) (Figure 10). According to Fujita et al. (2009), SnRK2 kinases are the main positive regulators of ABA signaling in response to water deficit stress in Arabidopsis.
The enriched terms exclusive to each treatment, ' AT' ('auxin transport') and 'PH' ('potassium ion homeostasis'), stand out for T25 (Figure 9). Regarding ' AT' , the scientific literature shows that auxin plays an essential role during abiotic stress-induced changes in the root. Developmental modifications to root system architecture (RSA) is vital for tolerance to drought and high salt stresses. It has long been established that RSA requires the phytohormone auxin. Auxin transport in plants also has received much attention [for a review, see Korver et al. (2018)]. The enrichment of the term 'PH, ' in turn, suggest that up-regulated VuPK may act in the referred process. Maintaining adequate plant potassium is critical for plant drought tolerance. A close relationship between these two factors has been demonstrated (Wang et al., 2013).
Finally, terms associated with brassinosteroid (a class of plant steroidal hormones), 'DBS' and 'BS, ' were specifically enriched from T150 up-regulated VuPKs (Figure 10). From the past decade, additional roles for brassinosteroids have emerged, particularly relating to stress tolerance (Ahammed et al., 2020). However, the two most modulated VuPKs for those terms, Q39011 (shaggy-related protein kinase) and Q9SCZ4 (receptor-like protein kinase FERONIA), are negative regulators of brassinosteroid response and brassinosteroid signal transduction pathway, respectively. Such observations suggest that this steroidal hormone does not actively participate in cowpea's response to root dehydration stress. However, further investigations to confirm this observation are necessary.
Our results provide a scientific foundation for the comprehensive understanding of VuPK omics aspects. In this work, we connect the VuPKs with the state-of-theart knowledge offering a general picture of this genomic presence, evolutionary aspects, and its participation in the cowpea response to biotic and abiotic stresses. Extensive works of functional characterization (transgenic or gene silencing approaches) will be the next steps to be performed to determine the specific roles (activation of signaling process or metabolic pathways, etc.) of VuPKs in unfavorable conditions. Comprehensive knowledge of these actors could aid the current knowledge about the cowpea molecular physiology. Additionally, these proteins could become a useful resource in the quest to generate stress-tolerant/resistant, high-yielding cowpea varieties.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https:// www.ncbi.nlm.nih.gov/, Root dehydration | BioProject ID: PRJNA605156; Mechanical injury + CABMV | BioProject ID: PRJNA655993; Mechanical injury + MV | BioProject ID: PRJNA656211.

AUTHOR CONTRIBUTIONS
AB-I conceived and designed the experiments. JF-N and AB performed the experiments. JF-N, MS, EK, and AB performed the laboratory experiments and analyzed the data. DM, GB, and JB-N provided support for the RNA-Seq assembly and quality analyses. AB-I obtained the financial resources. AB-I and JF-N wrote the manuscript. All authors read and approved the manuscript.

ACKNOWLEDGMENTS
The authors acknowledge the support of the Fundação de Amparo a Ciência e Tecnologia de Pernambuco (FACEPE), Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), and Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES) through the research fellowship grants.
Supplementary Table 1 | Reference genes and target VuPKs used in the relative expression analysis, presenting the primer pair sequences, melting temperature of the amplicon, and primer pair efficiency values.
Supplementary Table 2 | iTAK software output, indicating the VuPKs searched in the cowpea genome as well as the types of domain and their quantity per analyzed sequence, in addition to alignment values for the implemented Hidden Markov Model.
Supplementary Table 3 | Bona fide VuPKs with their respective kinase families as well as indication of the unknown kinases, whose identity was not corroborated by the two adopted categorization strategies.
Supplementary Table 4 | Characterization of the 118 kinase families of the cowpea genome, indicating the number of coding loci and 12 parameters that define the average gene structure of their bona fide VuPKs constituent.