Unraveling the Genomic Landscape and Stress Responses of Cowpea Thaumatin-Like Proteins: A Comprehensive Analysis

1. Introduction

In the course of the evolution, plants have developed a series of components to fight pathogens. Among several protein families involved in plant defense, we stand out thaumatin-like proteins (TLPs). TLPs are members of the Pathogenesis-related 5 (PR-5) gene family. This protein group received this name because of their structural similarity to proteins called thaumatins, isolated from the fruit of Thaumatococcus daniellii Benth. (Maranthaceae family) (van der Wel and Loeve, 1972). The mentioned group share the Thaumatin conserved domain (Pfam: PF00314) that covers almost 95% of the entire mature peptide (Petre et al. 2011). About their potential biological functions, TLPs have shown antifungal activity (Li et al. 2023; Wang et al. 2023; Yang et al. 2024; Liu et al. 2021), acting on fungal membrane permeabilization (Batalia et al. 1996). There is evidence that some TLPs representatives act on invading fungi through their binding activities or hydrolysis of β-1,3-glucans (Grenier et al. 1999; Trudel et al. 1998) or, still, through inhibition of enzymes such as xylanases (Fierens et al. 2007).

Hormones associated with plant immunity can also induce TLPs-encoding genes. For example, in tobacco seedlings, the induction of TLPs genes occurred in response to the application of ethylene (ET), methyl jasmonate (methylJA) and salicylic acid (SA) (Xu et al. 1994). A similar fact was also reported for transgenic rice, where the rice blast fungus elicitor, besides the hormones SA and methylJA, induced the β-glucuronidase (GUS) reporter gene driven by the Rtlp1 (Rice thaumatin-like protein 1) promoter (pRtlp1GUS) (Hiroyuki et al. 2008). TLPs expression are also associated with abiotic responses, as drought (Zhao et al. 2024; Muoki et al. 2021; Jung et al. 2005), wounding (Ruperti et al. 2002), and freezing (Liu et al. 2023; Hon et al. 1995). There are, additionally, reports that TLPs were induced after infection by virus or bacteria (Breiteneder and Radauer, 2004); however, in these situations, their action mechanism remains indeterminate. Regarding viruses, the work of Kim et al. (2005) observed that a tobacco (Nicotiana tabacum), TLP1, interacted with Cucumber mosaic virus proteins.

Due to their involvement in the response to pathogens and biotechnological potential, TLPs have been analyzed in a large number of species, such as Arabidopsis thaliana, rice (Oryza sativa), Populus spp., maize (Zea mays), Physcomitrella patens, Chlamydomonas spp., and wheat (Triticum aestivum) (Muoki et al. 2021; Sharma et al. 2022; Cao et al. 2016; Wang et al. 2010), among others. About legumes, information on this protein class is still quite limited (Feng et al. 2024). The vast majority of these studies, however, does not explicitly focus on global analysis of TLPs; only quote the some TLPs participation in the global transcriptional response, or focus on point TLPs. The only report found in the literature addressing the structural and functional genomics of TLPs was performed by Petre et al. (2011). These authors identified 42 validated TLP-loci in the genome Populus trichocarpa 'Nisqually-1' (version 2.0). In addition, it was observed that most TLPs were responsive to abiotic stresses (high ozone content, UV-B rays, drought, and high copper content) besides being also responsive to biotic stress (infection by the fungus Melampsora laricis-populina).

Cowpea is an important legume. Its social role is important: such crop is a relevant source of proteins and minerals for millions of people in sub-Saharan Africa and other developing regions (Boukar et al., 2016). Cowpea demonstrates exceptional phenotypic plasticity and strong physiological performance under a wide array of biotic and abiotic stresses. Given its importance and adaptive capacity, considerable efforts in omics research have been focused on this species. Recently, its reference genome was released (Lonardi et al., 2019), providing a vital resource for elucidating the crop’s physiological advantages under stress conditions, identifying genes with biotechnological potential, and serving as a key platform for comparative genomics within legumes. In Brazil, which contributes 12% of global cowpea production, the Cowpea Genomics Consortium (CpGC) was established. The CpGC manages sequenced genomes from multiple accessions and cultivars of this species, which are in various stages of assembly. In addition, it holds comprehensive transcriptomes (RNA-Seq) from cultivars demonstrating tolerance or resistance, which have been subjected to root dehydration or combined stressors (mechanical injury followed by viral inoculation with either Cowpea Aphid-borne Mosaic Virus or Cowpea Mosaic Severe Virus).

Thus, given the importance, multifunctionality, the responsiveness of plant TLPs under unfavorable conditions and the availability of CpGC data, the present work brings structural characterization and gene expression data of cowpea TLPs expressed under different stress types (mechanical injury, CABMV or CPSMV inoculation and root dehydration). The present work uncovered relevant information to the molecular physiology of cowpea under unfavorable conditions and provide data about the performance of TLPs in a legume, an issue still little addressed.

2. Materials and Methods

2.1. Biological Material, Experimental Design and Stress Application

❖

Root dehydration assay

Seeds of Vigna ungiculata cv. Pingo de Ouro (considered tolerant to water deficit and drought; Rodrigues et al., 2017; Bastos et al., 2011) were treated with 0.05 % (w/v) Thiram (tetramethylthiuram disulfide) and germinated during 2 days at 25° C ± 1 °C and 65% ± 5% of temperature and relative humidity, respectively. 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) in a randomized block experimental design, with three biological replicates (Supplementary Figure 1B). Each biological replicate was composed of two individuals. Plantlets were placed in supports in such a way that the roots of the seedlings were completely immersed in the solution (Supplementary Figure 1A). Plantlets were grown for three weeks (estágio de desenvolvimento V3) in greenhouse under a natural photoperiod of approximately 13/11 h light/dark cycle, temperature of 30 ± 5 °C and 60 ± 10 % relative humidity (RH). After this period, root dehydration treatment was initiated by withdrawing the nutrient solution from treated plants (Supplementary Figure 1A). Roots were collected after 25 min (RD25) and 150 min (RD150) after solution removing (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 1BB) were maintained in the nutrient solution and subsequently collected.

Throughout this manuscript, the expression contrast for ‘RD25 vs Cont.25’ was referred to as T25; similarly, the contrast for ‘RD150 vs Cont.150’ was referred to as T150.

❖

Mechanical injury and virus inoculation assays

The experiments involving mechanical injury and inoculation by CABMV (Cowpea Aphid-borne Mosaic Virus) or CPSMV (Cowpea Severe Mosaic Virus) were conducted under controlled conditions in a greenhouse at the Instituto Agronômico de Pernambuco (IPA; Recife, Pernambuco, Brazil) (Supplementary Figure 2A). The CABMV assay utilized the resistant genotype IT85F-2687 (Rocha et al., 1996; Oliveira et al., 2011), while the CPSMV assay employed the resistant genotype BR-14 Mulato (Cardoso et al., 1990).

The experimental procedures for both assays were conducted separately but followed similar protocols. Both accessions were sown and grown for three weeks, reaching the V3 development stage, under natural photoperiods and temperatures ranging from 28 to 32°C (Supplementary Figure 2A). The leaves of the youngest trifoliate were mechanically injured using carborundum (silicon carbide) to facilitate viral entry into the plant tissue, followed by the application of the viral inoculum (Supplementary Figure 2A).

Two post-injury/inoculation time points were implemented for each assay (Supplementary Figure 2B): 60 minutes and 16 hours. Each treatment had its corresponding absolute control or mock (Supplementary Figure 2B). The tissues were immediately frozen in liquid nitrogen and stored at -80°C until RNA extraction.

The experimental design was factorial, based on cultivar and post-inoculation period, with three biological replicates (BRs) for each control and treatment (Supplementary Figure 2B). Each BR consisted of five plants. All treatments were carried out in isolated areas to avoid cross-contamination through volatile compounds released by plants.

Differential gene expression in the two virus-related assays was induced by a combination of stress factors: mechanical injury and viral inoculation. Plant viruses cannot initiate infection independently - they require a vector or specific agricultural practices to breach the plant cell wall. According to Barna and Király (2004), plant viruses lack specific cellular receptors, unlike bacteriophages and animal viruses. Therefore, the combination of ‘mechanical injury and viral inoculation’ was employed to mimic natural infection processes.

Throughout this study, the expression contrasts for the mechanical injury (MI) and viral inoculation (CABMV or CPSMV) assays were defined as follows:

• MI+CABMV(60min) vs. control: MI_CABMV60’
• MI+CABMV(16h) vs. control: MI_CABMV16h
• MI+CPSMV(60min) vs. control: MI_CPSMV60’
• MI+CPSMV(16h) vs. control: MI_CPSMV16h

2.2. RNA Extraction and cDNA Synthesis

Total RNA was isolated using the “SV Total RNA Isolation System” kit (Promega, US) following manufacturer's protocol. Agarose gel (1.5 %) and Agilent 2100 Bioanalyzer (Agilent Technologies, EUA) evaluated both concentration and quality of total extracted RNA. Only samples with RNA integrity number (RIN) ≥ 8.0 were sequenced.

A “RNAm TruSeq^® Stranded LT-Set A” kit (RS-122-2101) (Illumina, San Diego, CA, EUA) 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 2500 system, using the following kits: “HiSeq^® Rapid PE Cluster Kit v2” (PE-402-4002); “SBS Kit v2” (200 Cycle; FC-402-4021); and “TruSeq^® Stranded mRNA LT - Set A” (RS-122-2101). All sequencing steps were performed at the Center for Functional Genomics, University of São Paulo (São Paulo, Brazil).

2.3. RNA-Seq Libraries Assembly and Differential Expression Analysis

The 12 RNA-Seq libraries sequenced for the root dehydration assay were combined with the 12 libraries from the “mechanical injury and viral inoculation by Cowpea Aphid-borne Mosaic Virus (CABMV)” experiment and the 12 libraries from the “mechanical injury and inoculation by Cowpea Severe Mosaic Virus (CPSMV)” experiment. This combined assembly facilitated the generation of longer and more robust transcripts. A total of 72 replicates (12 biological replicates for root dehydration, 12 for CABMV, and 12 for CPSMV, each with two technical replicates) were processed through the pipeline, briefly outlined below.

Initially, the Trimmomatic v0.36 package was used to remove adapters and low-quality sequences (Bolger et al., 2014). Reads were trimmed from the 3' end to ensure a Phred score of at least 30. Illumina sequencing adapters were removed, and only reads with a minimum length of 32 bp were retained. Subsequently, de novo transcriptome assembly was performed using the Trinity software (Grabherr et al., 2011), as described by Haas et al. (2013). The quality of the assembled transcriptomes was assessed by evaluating the N50 of the resulting libraries, along with the count and distribution of contig lengths.

Notably, differential expression analysis was conducted independently for each assay (root dehydration; mechanical injury + CABMV inoculation; mechanical injury + CPSMV inoculation). This analysis was performed using the edgeR tool (Robinson et al., 2010) implemented within the Bioconductor package (Huber et al., 2015). Transcripts with Log2FC values between -1 and 1, p < 0.05, and FDR < 0.05 were considered differentially expressed.

Hierarchical clustering of VuTLPs was carried out with CLUSTER 3.0 (Eisen et al. 1998), and the heatmap was visualized using TreeView (Page, 1996). In addition, Venn diagrams were also generated using software Venny (Oliveros, 2007).

2.4. Cowpea TLPs (VuTLPs) in Silico Characterization

The in silico analysis was carried out using reads of cowpea RNA-Seq libraries under CABMV or CPSMV inoculation or under root dehydration (Supplementary Table S1). VuTLPs were predicted based on sequence homology searches, using thaumatin-domain with sequences previously characterized (see Jesús-Pires et al. 2019). An outline of the annotation strategy is presented in Fig. 1. For the identification of VuTLPs candidates, the alignments were carried out against CpGC database using tBLASTn (Altschul et al. 1990) (Figure 1) with a cut-off value of 1e⁻⁰⁵. VuTLP candidates were annotated against NCBI and UniProt and analyzed for their score, e-values, sequence size, and presence of Conserved Thaumatin Domains (CTDs), as shown in Supplementary Table S2.

All sequences were screened for conserved motifs, with the aid of rps-BLAST on CD-search tool (Altschul et al. 1990) (Figure 1). Only orthologs presenting the expected features (CTDs and motifs) were translated into proteins and considered for subsequent evaluation (Figure 1).

Aiming to enrich the characterization of predicted VuTLPs, the identified amino acid sequences were further characterized including the following analyzes (Figure 1):

• Determination of the putative isoelectric point and molecular weight using ProtParam (https://web.expasy.org/protparam/; (Gasteiger et al. 2005); Figure 1);
• Prediction of signal peptide for each cowpea TLP candidate with SignalP 4.1 Server (Petersen et al. 2011; http://www.cbs.dtu.dk/services/ SignalP/);
• Prediction of the subcellular localization with TargetP 1.1 Server (Emanuelsson et al. 2007; http://www.cbs.dtu.dk/services/TargetP).

Additionally, we used Clustal Omega (http://www.ebi.ac.uk/Tools/msa/clustalo/; Sievers et al. 2011) to generate a multiple sequence alignment of full-length VuTLPs compared with previously described TLPs from other plant species, aiming at the localization and comparison of CTDs and residues among TLP sequences. For this step, some homolog proteins denominated “osmotin” and “zeamatin” were also employed (Supplementary Figure 3; Supplementary Material SM1), considering literature evidences (Singh et al. 1989; Roberts and Selitrennikoff, 1990; Skadsen et al. 2000) that such proteins present similarities with TLPs.

2.5. Distribution of VuTLP-Candidates in V. unguiculata Genome

To determine the genomic distribution of VuTLPs, cowpea candidates were aligned against V. unguiculata genome V12 available at Phytozome Database (https://phytozome.jgi.doe.gov/) aiming to map these sequences in virtual chromosomes through the BLASTp tool. This step aimed to infer on the distribution, relative position, and abundance of TLPs-coding loci. We considered the best-hit (e-value cut off < 1e⁻¹⁰) to allow the identification of VuTLPs along the virtual chromosomes. Afterward, the identified mapping positions were used to build a virtual ideogram (Figure 2) regarding the TLPs distribution on V. unguiculata pseudochromosomes (n = 11).

2.6. NJ Analysis and VuTLPs Gene Features

For the Neighbor-Joining (NJ) analysis, used 34 TLP-coding genes. The multiple alignment was submitted to manual edition and elimination of not aligned extremities and uninformative autapomorphies. The NJ analysis was carried out using the MEGA 7 (Kumar et al. 2016). Statistical support of the branches was evaluated by a bootstrap analysis with 2,000 replicates.

The gene structure was examined in terms of exon-intron organization of each VuTLP. The display of exon-intron was performed using GSDS 2.0 server (Hu et al. 2015). The VuTLPs gene features também foram determinadas via Genestats script.

The following parameters were scrutinized:

1. transcript sequence length; 2. number of exons; 3. total exon sequence length; 4. number of introns; 5. total intron sequence length; 6. number of CDS chunks; 7. total CDS sequence length; 8. number of 5' UTR sequences; 9. total 5' UTR sequence length; 10. number of 3' UTR sequences; 11. total 3' UTR sequence length.

2.7. Analysis of VuTLPs Duplication

Multiple Collinearity Scan toolkit (MCScanX; Wang et al., 2012) package - downstream analysis mode: ‘duplicate_gene_classifier’ - was applied to classify the origins of the duplicate VuTLPs genes of the cowpea genome, into:

• Whole genome /segmental (i.e., collinear genes in collinear blocks);
• Tandem (adjacent loci in a genome region);
• Proximal (gene pairs in nearby chromosomal region but not adjacent);
• Dispersed (other modes than segmental, tandem and proximal) duplications.

The following procedure was used by MCScanX to assign the duplication mechanisms:

a)

All genes were initially classified as ‘singletons’ (i.e., without duplicates in the genome) and assigned gene ranks according to their order of appearance along chromosomes;

b)

BLASTp results were evaluated and the genes with BLASTp hits to other genes were re-labeled as ‘dispersed duplicates’;

c)

In any BLASTp hit, the two genes were re-labeled as ‘proximal duplicates’ if they had a difference of gene rank < 20 (configurable);

d)

In any BLASTp hit, the two genes were re-labeled as ‘tandem duplicates’ if they had a difference of gene rank = 1;

e)

MCScanX was executed. The anchor genes in collinear blocks were re-labeled as ‘WGD/segmental’.

f)

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.

2.8. Ratio of Synonymous (Ks) and Non-Synonymous (Ka) Substitutions per Site for Tandem Duplicated Genes

The ratio of Ka to Ks (Ka/Ks) was used to determine the selection pressure among tandem and segmental duplication events. ClustalW 2.0 (Larkin et al., 2007) software first aligned the full-length coding sequences of tandemly and segmentaly duplicated VuTLPs genes. Then, the Ka and Ks 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 Ka/Ks ratio is an indicator of the selection history on genes or gene regions:

• if its value is lower than ‘1’, the duplicated gene pairs may have evolved from purifying selection (also called as negative selection; conserves the amino acid sequence);
• If Ka/Ks equal ‘1’, means neutral selection (that had no constraint for sequence divergence);
• while Ka/Ks greater than ‘1’ means positive selection (that led to different peptides).

2.9. Promoter Analysis

Cowpea promoter regions (up to 1.0 kb upstream of TLP-coding genes) were downloaded from Phytozome database v8.0 (https://phytozome.jgi.doe.gov). The motifs (candidate cis-regulatory elements) in each promoter were revealed through the MEME v5.0.3 software (Bailey and Elkan, 1994) (http://meme-suite.org/tools/meme). The software reports an e-value for each identified motif and gives an estimate of the number of expected motifs found by chance if the input sequences were shuffled. In the present simulation, an e-value < 10^-2 was adopted as the cut-off for characterization of bona fide cis-regulatory element candidates. The maximum number of motifs searched in the present analysis was 10. The motifs analyzed exhibited between six and 50 nucleotides in length.

After MEME software analysis, the Tomtom v4.11.2 software (Gupta et al. 2007) (http://meme-suite.org/tools/tomtom) was used with the JASPAR database (file: JASPAR2018_CORE_plants_non-redundant) to annotate putative cis-regulatory elements. JASPAR searches one or more queries (candidate cis-regulatory element) against annotated motifs ranked by p-value (cut-off < 10^-2). The q-value (false discovery rate; cut-off < 10^-2) of each match is also reported. In this work, the presented identities of the cis-regulatory element candidates were associated with the target motif that exhibited the most significant p-value.

2.10. qPCR Analysis

For validation of the RNA-Seq data of the different experiments, RT-qPCR analyses were performed with treatments and control for each time sampled. For each time point three biological replicates, each with three technical replicates, were analyzed to ensure statistical reliability. The qPCR reactions were performed on the CFX96 Touch Real-Time PCR Detection System (Bio-Rad), using SYBR Green detection. The PCR program was adjusted to 95 °C for 2 min, followed by 35 cycles of 95 °C for 15 s, 60 °C for 1 min, and 72 °C for 15 s. After amplification, dissociation curves were produced (65–95 °C at a heating rate of 0.5 °C / sec and acquiring fluorescence data every 0.3 °C) to confirm the specificity of the PCR products.

The amplification efficiency (E = 10^{(−1/slope of the standard curve)} for all primer pairs was

determined from a 4-point standard curve generated by serial dilutions of cDNA (10-fold each) in technical triplicates. β-tubulin-2, F-box and UBQ10 were used as reference genes for normalization of MI_CABMV and MI_CPSMV assays while Actin and UE21D were used as reference genes for normalization of radicular dehydration assay (Amorim et al. 2018). Primers for VuTLPs candidates were designed based on the V. unguiculata transcriptome using the Primer3 tool (http://bioinfo.ut.ee/primer3-0.4.0/) under program´s default settings (Supplementary Table S3).

The Rest2009 (Pfaffl et al. 2009) software package (standard mode) was used for relative expression analysis of target transcripts. Hypothesis testing (p < 0.05) was used to determine whether the differences in target transcripts expression between the control and treatment conditions were significant.

3. Results

3.1. In Silico Identification and Characterization of VuTLPs

The tBLASTn search allowed the identification of 126 putative VuTLPs (Supplementary Table S2) expressed under biotic and abiotic stresses. From these 33 presented no ORF associated with TLP, while 32 presented complete TLP conserved domain (CTD) and 61 presented incomplete CTDs. Among candidate VuTLPs with complete CTD, 30 were expressed in cowpea RNA-Seq libraries generated under radicular dehydration, while 28 and 27 VuTLPs were expressed under CABMV and CPSMV inoculation, respectively. Considering the total amount (abiotic + biotic stresses), 30 non-redundant VuTLPs were recognized (Supplementary Material SM2).

All 30 predicted VuTLPs presented CTD in at least one searched database (NCBI CDD: cd09218 and cl02511, Pfam: pfam00314, SMART ID: smart00205; Supplementary Table S2). Five highly conserved amino acids (REDDD motif), in an acidic cleft, were detected in most VuTLPs. However, the alignment of the referred candidates with other TLPs from A. thaliana, Nicotiana tabacum, Zea mays, and Cicer arietinum revealed that some VuTLPs did not contain all expected amino acids in REDDD motif. For this group, some VuTLPs exhibited the following substitutions (Supplementary Figure 3): the glutamic acid (E) by a glutamine (Q) in the transcript Vu54351|c0_g1_i1 and in the transcripts Vu85277|c0_g1_i1 and Vu28453|c0_g1_i1 the glutamic acid (E) was replaced by aspartic acid (D) of the same functional group; and the first aspartic acid (D) by a glycine (G) in the transcript Vu54351|c0_g1_i1. Other modifications in the mentioned motif were also perceived. For instance, in the transcripts Vu85277|c0_g1_i1 and Vu28453|c0_g1_i1 (Supplementary Figure 3) the second aspartic acid (D) was replaced by an asparagine (N).

Regarding the 16-cysteine conservation, just three cowpea sequences (Vu58024|c0_g1_i1, Vu49480|c1_g2_i2 and Vu49480|c1_g2_i1; Supplementary Figure 3) lacked the first cysteine; whereas the sixth cysteine of transcript Vu83475|c0_g1_i1 was replaced by a serine (S). In turn, all other VuTLPs presented the 16 residues. Additionally, expressed VuTLPs exhibited a molecular weight between 18.97 and 34.93 kDa and an isoelectric point varying from 4.23 to 9.05 (Supplementary Table S2).

The SignalP prediction, in turn, revealed that 20 VuTLPs present hydrophobic signal peptide sequences.

3.2. Mapping of VuTLPs in V. unguiculata Genome and Analysis of Gene Duplication

The approach revealed the distribution of the VuTLPs in 34 loci located along the chromosomes, with exception chromosome ten (Figure 2).

A clustering of VuTLPs loci was observed, being more pronounced on chromosomes two and eight (Figure 2; Supplementary Table S4). For example, the long arm of V. unguiculata chromosome 2 presents six and two clustered TLP coding genes (Figure 2). Such a distribution is indicative of evolutionary processes through tandem duplication. Besides the observed clustering in most chromosomes, the position of VuTLP coding genes was mainly subterminal or intercalary, while a pericentric position of VuTLPs is observed only in chromosomes 6 and 7 (Figure 2).

Gene duplications, including segmental and tandem duplication, have been considered as one of the main forces in the evolution and expansion of a gene Family (Lynch and Conery, 2000). The analysis of the expansion of TLP-coding genes revealed that the 34 TLP-coding genes have expanded through gene duplication processes (Supplementary Table S5). In this study, we found four types of gene duplications: dispersed, proximal, tandem, and WGD or segmental duplication by MCScanX program. Our analysis provides evidence that whole-genome duplication (WGD) or segmental duplication has been the primary mechanism driving the expansion of VuTLPs in the cowpea genome, resulting in a total of 15 gene pairs (Supplementary Table S6). Moreover, VuTLP genes were observed to cluster within six tandem duplication regions across four of the 11 chromosomes (Supplementary Table S7). Finally, we identified five dispersed duplication events and one proximal duplication event contributing to the expansion of these genes (Supplementary Table S5).

3.3. Ratio of Synonymous and Non-Synonymous Substitutions for Tandem Duplicated Genes

The forces that drive natural selection can be inferred from the types of nucleotide substitutions in coding sequence of the genes. An informative parameter of the evolution of genes under selection has been the ratio of the rate of nonsynonymous substitution (Ka; that causes an amino acid change) to the rate of synonymous substitution (Ks; that not causes an amino acid change). Os valores de Ka/Ks variaram de 0,00 a 1,36 for tandem duplication (Supplementary Table S7) e de 0,00 a 0,31 for segmental duplication (Supplementary Table S6). Sixteen pairs (16 duets) presented a Ka/Ks ratio of less than 1 (Supplementary Tables S6-S7), suggesting that these genes have been subjected to purifying selection. Em outro contexto, observou-se que cinco duetos apresentaram Ka/Ks > 1 (Supplementary Table S6-S7), implying that positive selection played an important role in the evolution desses duetos.

3.4. NJ Analysis and Exon-Intron Organization

A Neighbor-Joining tree was constructed using the protein sequences of all the 34 VuTLP genes. The cowpea VuTLP family was divided into three groups (Supplementary Figure 4). The Group III colored in purple (58.82%) contained the most members, followed by Group II colored in pink (27.27%) and the least represented group was the I colored in green (5.88%) with only two VuTLPs (VuTLP16 and VuTLP32). To understand the formation of gene groups, several characteristics were analyzed, including the presence of a signal peptide, molecular weight, isoelectric point, presence of REDDD residues, and subcellular localization (Table 1). Group I comprised genes with molecular weights ranging from 24.16 to 25.18 kDa and isoelectric points from 6.80 to 7.32. Notably, in the gene VuTLP32, alterations were observed in the REDDD residues, where E was replaced by Q, and the first D by G. Group II included genes with molecular weights ranging from 21.48 to 29.20 kDa and isoelectric points from 4.04 to 8.76. In the gene VuTLP12, a change in the REDDD residues was noted, with the second D replaced by N. Group III consisted of genes with molecular weights from 13.75 to 108.18 kDa and isoelectric points from 4.00 to 8.19. In the gene VuTLP6, the REDDD residue showed an alteration where R was replaced by Q. Additionally, in genes VuTLP7, VuTLP5, and VuTLP18, the first D residue was substituted by Y, and in VuTLP18, the second D was replaced by Q. Moreover, in this group, two genes (VuTLP5 and VuTLP18) were localized to the plasma membrane.

Besides, the most of thirty-four VuTLP genes shown in gene structure analysis had three exons, while VuTLP7, VuTLP31, VuTLP20, VuTLP21, VuTLP29, VuTLP1, VuTLP17, VuTLP15, VuTLP24, VuTLP25 and VuTLP16 had two and VuTLP4, VuTLP19, VuTLP9, VuTLP11, VuTLP8, VuTLP10 and VuTLP32 had one exon. Four genes (VuTLP11, VuTLP8, VuTLP10 and VuTLP32) had no introns (Supplementary Figure 4).

3.5. Promoter VuTLPs

The 34 cowpea TLP-coding genes had their predicted promoters (1 kb) analyzed regarding the presence and identity of candidate cis-regulatory elements (CCREs), via the MEME and Tomtom tools. Four of the identified promoter motifs were within the stipulated cut-off (e-value < 10^-2, Figure 3; Supplementary Figure 5A, 5B, 5C e 5D), being considered as bona fide CCREs, (see Supplementary Figure 5A, 5B, 5C e 5D). The six remaining motifs (Figure 3) were associated to TFs [Dof-type (dark green and pink boxes), C2H2 (light green), G2-like (orange), and CPP (light blue), and one without specific annotation] but did not achieve the specified parameters, being excluded from the discussion.

The identified bona fide CCREs are associated with the three TF families (Figure 3): bHLH (two motifs; light green and blue boxes, Supplementary Figure 5A and 5C), MYB-related (light blue, Figure 3; Supplementary Figure 5B) and Dof-type (red, Figure 3; Supplementary Figure 5D). CEMR associated with a Dof-type transcription factor (JASPAR ID MA1267.1; Figure 3; Supplementary Figure 5D) was the most abundant promoter, observed in 42 sites, followed by MYB-related (JASPAR ID MA1182.1; Figure 3; Supplementary Figure 5B) and bHLH (JASPAR ID MA0587.1; Figure 3; Supplementary Figure 5C) each anchored in 11 sites. Finally, bHLH (JASPAR ID 1361.1; Figure 3; Supplementary Figure 5A) anchored in five sites.

3.6. TLP Content and Expression in Cowpea Transcriptomes under Different Stress Types

Considering all available RNA-Seq libraries for biotic and abiotic stresses (Supplementary Table S8), 30 VuTLPs were identified. A partir daqui iremos dividir os ensaios para melhor analisar a regulação das VuTLPs.

3.6.1. Mechanical Injury (MI) and Virus Inoculation Assays

Twenty-eight VuTLPs-coding transcripts were expressed in the MI_CABMV RNA-Seq libraries (Supplementary Table S8; Supplementary Figure 6A). Of these, six (Vu3908|c0_g1_i1, Vu93958|c3_g1_i3, Vu85277|c0_g1_i1, Vu95346|c0_g1_i1, Vu58024|c0_g1_i1 and Vu83475|c0_g1_i1) were up-regulated exclusively at 60 min (Supplementary Figure 6B), indicating that VuTLPs are preferably recruited during the first hours after the MI_CABMV treatments. In addition, only one VuTLP-encoding transcript (Vu44758|c0_g1_i1) showed up-regulation at 16 h. Besides, the transcript Vu9781|c2_g1_i3 encoding VuTLP was up-regulated at 60 min and repressed at 16 h treatment, whereas the others presented constitutive expression or were not expressed in the comparisons.

Regarding the MI_CPSMV assay, twenty-seven VuTLPs-encoding transcripts were expressed (Supplementary Table S8; Supplementary Figure 7A). Of these, five were up-regulated exclusively at 60 min (Vu95899|c0_g1_i1, Vu133945|c0_g1_i1, Vu9781|c2_g1_i3, Vu9763|c0_g1_i1 and Vu95346|c0_g1_i1) and three at 16 h (Vu49480|c1_g2_i2, Vu88791|c0_g1_i2 and Vu44758|c0_g1_i1) (Supplementary Figure 7B). In addition, a VuTLP-encoding (Vu83475|c0_g1_i1) transcript was up-regulated in both treatments (60 min and 16 h; Supplementary Figure 7A). These results also point out that potential VuTLPs were recruited in the first hour after MI_CPSMV treatment.

When comparing the transcriptional orchestration in response to the two MI_viral inoculation assays, it was observed that four potential VuTLPs (Vu9781|c2_g1_i3, Vu95346|c0_g1_i1, Vu44758|c0_g1_i1 and Vu83475|c0_g1_i1) were up-regulated in response to both assays (Supplementary Figure 8A). Por sua vez, four (Vu3908|c0_g1_i1, Vu93958|c3_g1_i3, Vu85277|c0_g1_i1 and Vu58024|c0_g1_i1) and five (Vu95899|c0_g1_i1, Vu133945|c0_g1_i1, Vu9763|c0_g1_i1, Vu49480|c1_g2_i2 and Vu88791|c0_g1_i2) VuTLPs-encoding transcripts were induced exclusively in response to MI_CABMV and MI_CPSMV, respectively (Supplementary Figure 8A).

3.6.2. Root Dehydration Assay

Thirty potential VuTLPs-encoding transcripts were expressed in cowpea under root dehydration (Supplementary Table S8; Supplementary Figure 9A). Of these, Vu73314|c0_g1_i1 and Vu49480|c1_g2_i1 were up- and repressed at 25 min, respectively. Additionally, the transcript Vu94809|c0_g2_i1 encoding VuTLP was induced at 25 min and repressed at 150 min treatment. In addition, the transcript Vu88791|c0_g1_i2 was induced in both treatments (25 min and 150 min). Besides, eight (Vu88791|c0_g1_i1, Vu8524|c1_g3_i1, Vu93958|c3_g1_i3, Vu93958|c3_g1_i2, Vu49480|c1_g3_i1, Vu9781|c2_g1_i3, Vu149510|c0_g2_i1 and Vu83475|c0_g1_i1) and three (Vu54351|c0_g1_i1, Vu3908|c0_g1_i1 and Vu79818|c1_g2_i8) transcripts encoding for VuTLPs were induced and repressed exclusively at 150 min treatment, respectively. These results indicate that some TLPs are recruited in the first minutes after stress imposition, even though a higher number of transcripts were detected at 150 min (Supplementary Figure 9B).

3.6.3. MI_CABMV vs. MI_CPSMV vs. Root Dehydration Assays

A comparison between the transcriptional orchestration in response to both MI_viruses (in leaves) and root dehydration (roots) revealed interesting insights. There were 20 potential VuTLPs expressed jointly in the assays (MI_CABMV and MI_CPSMV and root dehydration, Supplementary Table S8; Supplementary Figure 8B), thus indicating a crosstalk role for TLPs under different stress types and tissues.

Qualitative analyses showed a heterogeneous orchestration of VuTLPs. Only two VuTLPs candidates (Vu9781|c2_g1_i3 and Vu83475|c0_g1_i1) were up-regulated in the three assays (Supplementary Table S8). Exclusive responses were also observed under different treatments. For example, three (Vu3908|c0_g1_i1, Vu85277|c0_g1_i1 and Vu58024|c0_g1_i1), four (Vu95899|c0_g1_i1, Vu133945|c0_g1_i1, Vu49480|c1_g2_i2 and Vu9763|c0_g1_i1) and seven (Vu88791|c0_g1_i1, Vu8524|c1_g3_i1, Vu94809|c0_g2_i1, Vu93958|c3_g1_i2, Vu49480|c1_g3_i1, Vu73314|c0_g1_i1 and Vu149510|c0_g2_i1) VuTLPs-encoding transcripts were up-regulated exclusively under MI_CABMV and MI_CPSMV inoculation and root dehydration, respectively. Thus, there is no overall conservation in VuTLPs transcriptional orchestration towards different analyzed assays.

3.7. RNA-Seq Data Validation by qPCR

3.7.1. MI_CABMV and MI_CPSMV Assays

The selection of target transcripts (Vu44758|c0_g1_i1, Vu9781|c2_g1_i3, Vu83475|c0_g1_i1, Vu93958|c3_g1_i3 and Vu58024|c0_g1_i1) was based on their regulation in the RNA-Seq libraries (Supplementary Table S8) for MI_CABMV experiment. Considering the functional primer pairs, todos apresentaram eficiência dentro dos padrões aceitáveis (Supplementary Figure 10). Relative expression analysis showed that Vu9781|c2_g1_i3, Vu83475|c0_g1_i1, Vu93958|c3_g1_i3 and Vu58024|c0_g1_i1 were up-regulated (UR) at 60 min (Supplementary Figure 11) while Vu44758|c0_g1_i1 was UR at 16 h (Supplementary Table S9 and Supplementary Figure 11), confirming RNA-Seq data. However, qPCR results have shown that Vu9781|c2_g1_i3 was not modulated at 16 h (not significant; ns), differing from RNA-Seq results where it was down-regulated (Supplementary Table S9 and Supplementary Figure 11).

For MI_CPSMV treatment, three differentially expressed target transcripts (Vu9781|c2_g1_i3, Vu88791|c0_g1_i2 and Vu44758|c0_g1_i1) were evaluated, todos os pares de primers analisados nesta etapa apresentaram eficiência dentro dos padrões aceitáveis; Supplementary Figure 12). Vu9781|c2_g1_i3 was UR at 60 min and Vu88791|c0_g1_i2 and Vu44758|c0_g1_i1 were UR at 16 h (Supplementary Table S9 and Supplementary Figure 13). All these results confirmed RNA-Seq expression data.

3.7.2. Root Dehydration

For the mentioned assay, six differentially expressed target transcripts (Vu88791|c0_g1_i2, Vu9781|c2_g1_i3, Vu3908|c0_g1_i1, Vu54351|c0_g1_i1, Vu79818|c1_g2_i8 and Vu49480|c1_g3_i1) were analyzed by qPCR, all with acceptable efficiency values (93.00 to 107.70%; Supplementary Figure 14). Relative expression analysis showed that Vu88791|c0_g1_i2, Vu9781|c2_g1_i3 and Vu49480|c1_g3_i1 were up-regulated at 150 min and transcripts Vu3908|c0_g1_i1 and Vu54351|c0_g1_i1 were down-regulated at 150 min (Supplementary Table S9; Supplementary Figure 15). All these results confirmed RNA-Seq expression data.

4. Discussion

Using TLPs sequences of others species as seed sequences, cowpea VuTLPs candidates were searched in the RNA-Seq libraries generated under radicular dehydration or CABMV or CPSMV inoculation. Such conditions were chosen to exemplify the participation of TLPs in response to both, biotic and abiotic stresses, since TLP participation in response to abiotic stress is still incipient when compared to its biotic counterpart. Additionally, RNA-Seq reads have up to 400 bp, depending on the sequencing technology used (Wang et al. 2009), allowing direct VuTLP structural characterization.

Most of the identified VuTLPs presented the REDDD motif, in an acidic cleft. This configuration is responsible for the reported antifungal activity of TLPs in plants (Cao et al. 2016; Ghosh and Chakrabarti, 2008). However, the alignment of the referred candidates with other TLPs from A. thaliana, N. tabacum, Z. mays, and C. arietinum revealed that some VuTLPs did not contain all expected amino acids in REDDD motif. It is still not clear whether these small differences have a significant impact on the substrate selectivity or protein function (Petre et al. 2011). Petre et al. (2011) found that the acidic cleft is the most conserved region among the eukaryotic TLPs. For example, in sTLPs (small TLPs), the most of REDDD amino acids are conserved, but wheat sTLPs sequences with resolved structure have revealed no acidic clefts nor any particular conserved region which may be linked to the xylanase inhibitor function reported for such proteins (Petre et al. 2011, Fierens et al. 2007).

Regarding the 16-cysteine conservation, with the exception of most VuTLPs, all other VuTLPs presented the 16 residues that form eight disulfide bonds necessary for correct folding and to ensure a high level of thermostability and pH constancy (Liu et al. 2012). Most TLPs presented molecular weights ranging from 21 to 26 kDa (Liu, Sturrock and Ekramoddoullah, 2010) and isoelectric point ranging between 3.4 and 12.0 (Velazhahan, Datta Muthukrishnan, 1999), having this last one influence on the electrostatic potential at the molecular surface (Tachi et al. 2009). Thus, our results show that most of the VuTLPs presented values within the expected standard for both characteristics, except for some VuTLPs that presented specificities regarding their molecular weight.

The SignalP prediction revealed that most VuTLPs present hydrophobic signal peptide sequences, indicating that they are predominantly secreted, a result confirmed by TargetP prediction. TLPs are located mainly in two subcellular compartments; apoplastic space or vacuole because of an N-terminal signal peptide that directs the mature protein to a secretory pathway and in some instances, a C-terminal polypeptide that directs them to a vacuolar compartment (Melcher et al. 1993). Previous studies have used GFP-gene approach associated with a transient expression system to study the localization of these proteins. Most of them point to a location in extracellular spaces (Singh et al. 2013; Zhang et al. 2017; Yang et al. 2024), corroborating with the in silico prediction.

To infer on the distribution, relative position, and abundance of VuTLPs, we anchored candidate transcripts against the cowpea genome, recently available (Lonardi et al. 2019). The amount of VuTLPs is similar to those found in other leguminous plant genomes. Specifically, Ammopiptanthus nanus possesses 31 TLP genes, Vitis vinifera has 32, Medicago truncatula has 49, Trifolium pratense has 23, and Lupinus albus has 29 (Liu et al., 2023).

The results of gene duplication analysis suggested the existence of gene replication events in cowpea genome. VuTLP genes were found to cluster into six tandem duplication event regions on 4 of 11 chromosomes, implying that they may originate from the recent gene duplication events. Moreover, segmental duplications were common in the cowpea TLP family. Fifteen segmental duplication gene pairs were identified in cowpea genome, compared to 3, 2, 7 and 4 in A. thaliana, rice, poplar and maize, suggesting that a majority of the VuTLP genes were generated by segmental replications (Cao et al. 2016). Cao et al. (2016) suggested that tandem and segmental duplications are the main factors leading to the enlargement and diversity of TLP family in plants.

Yang et al. (2024) identifed six segmentally duplicated TLP gene pairs in the Panax notoginseng genome. Gene duplication analysis revealed thirteen gene duplication events in the watermelon genome, including one tandem and twelve segmental duplications (Ram et al. 2022). Liu et al. (2020b) identifed twelve segmentally duplicated and seven tandemly duplicated gene pairs that have experienced purifying selection in the melon genome.

Natural selection analysis is widely used to estimate how natural pressure and evolutionary forces affect duplicated genes and corresponding proteins (Ellegren, 2008). In cowpea, a majority of TLP genes were shown to be driven by purifying selection, but few sites of TLP sequences were shown to be driven by positive selection. However, very similar results were found in the study on the evolution of TLPs in A. thaliana, rice and maize, which revealed that many TLP sequences were under purifying selection pressure and only on a few sites of TLP sequences may have experienced positive selection during the evolutionary process (Cao et al. 2016). Study on Poplar TLP genes suggested that four TLP kinases and ten TLPs have undergone some positive selection, indicating that plant TLPs may have experienced a diversified natural selection processes (Liu, Sturrock and Ekramoddoullah, 2010). Some of the exposed amino acids of TLPs are under positive selection, the amino acids that form the acidic cleft are always under purifying selection, which indicated that conservation of the acidic cleft are crucial and could be important for the antifungal activity of TLPs (Misas-Villamil and van der Hoorn, 2008).

In the current study, 34 cowpea VuTLP genes were analyzed into three groups through the NJ tree. In all obtained groups, a clear distinction emerged regarding the VuTLP genes when analyzing the REDDD motif associated with duplication events. Group III, which contained the highest number of clustered genes, displayed a consistent pattern of residue substitution in the REDDD motif, where the first D was replaced by Y. In contrast, genes in the other groups exhibited different substitution patterns. When comparing the conservation and substitution patterns of the REDDD residues with the duplication events, it was observed that residues maintaining conservation were expanded through segmental duplication and were under purifying selection. This indicates that these residues, and consequently the antifungal activity associated with this motif, were preserved. Conversely, residues that underwent alterations or substitutions expanded through tandem duplication and were subject to positive selection, suggesting a diversification of these residues that could lead to the emergence of new functions. Analyzing eukaryotic TLPs, Liu, Sturrock and Ekramoddoullah (2010) observed the formation of nine distinct groups showing, what he called TLP superfamily, as highly diversified. A possible explanation for the formation of these groups can be attributed to the structural diversity presented by TLPs. According to Petre et al. (2011), such diversity may influence the biological and biochemical functions, and differences in the topology around the cleft could determine the specificity of TLPs to their target ligands (Min et al. 2004). This may be true also considering the VuTLPs since their diversity and position in different positions of the NJ tree is possibly due to the presence of orthologs and paralogs, with structural and functional specificities. Zhao et al. (2024) identified four distinct groupings within the TLP gene Family and highlighted the presence of distinct subfamilies within each group indicating the close evolutionary associations between certain TLP genes from faba bean and other legume species.

Over time plants acquired new weapons to promote resistance over time (Burdon and Thrall 2009; Benko-Iseppon et al. 2010) while pathogens also continued to evolve, bringing pressure on the structure and diversity of some associated gene families, increasing the size of some protein families (Kohler et al. 2008), including PR-genes. Such an increase represents a major significance for the functional diversity via sub- or neo-functionalization of paralogs (Paterson et al. 2010). In this scenario, the natural selection of genes with new functions under environmental pressure probably played a significant role also in the evolution of TLPs (Liu, Sturrock and Ekramoddoullah, 2010), as may also be suggested for cowpea.

In addition, regarding the gene structure of cowpea TLPs, genes with three, two and one exon were verified, besides to four genes that did not present introns. The organization of introns and exons of TLPs have been reported to be in the range of one to ten exons (Liu et al., 2010; Petre et al., 2011; Cao et al., 2016). Out of the 10 TLP genes identified in the faba bean transcriptome, two genes exhibited three exons, five genes consisted of only one exon, and three genes displayed two exons (Zhao et al. 2024). Yan et al. (2017) found that most of the grape TLP genes had less than four exons and four genes also had no introns.

The 34 cowpea TLP coding genes had their predicted promoters (1 kb) analyzed in regard to the presence and identity of CCREs. Considering CCREs, a MEME motif is a sequence pattern that repeatedly occurs in one or more sequences in the input group. MEME can be used to discover novel patterns because discoveries are based only on the input sequences, not on any prior knowledge, such as databases of known motifs (Bailey and Elkan, 1994). Discovered bona fide CCREs can be used to identify potential physiological or adaptive processes in which TLPs participate, associating transcription factors and their intrinsic biological processes. This action adds value to the molecular dynamics of cowpea TLPs, especially in association with the considered transcriptomic data.

Concerning the promoters analyzed, CCREs associated with bHLH, Dof-type, and MYB-related TFs attended the stipulated statistical parameters. Also considering literature data, these TFs present an intimate association with plant stress response and other physiological processes. bHLH members can act as transcriptional activator or repressor and play essential roles in metabolic and developmental processes (Feller et al. 2011). Some bHLH have also been reported to be modulated under cold, drought, and salt stress (Wang et al. 2018; Mao et al. 2017). Dof-type transcription factors, in turn, participate widely in plant development and abiotic stress response (Li et al. 2016; Wang et al. 2017). In tobacco, the Sar8.2b gene can be activated by the Dof-type TF, which is related to systemic acquired resistance (Song and Goodman, 2002). Finally, there are reports on MYB-related TFs acting as retrotransposon regulators and controlling defense-related genes, besides being induced by wounding and other elicitors, in tobacco plants (Sugimoto, 2000). Such in silico data on promoters associated with TLPs in RNA-Seq libraries under different stresses (biotic and abiotic) indicate a plurality of roles for TLPs, also reinforcing their biotechnological potential.

The RNA-Seq data showed that potential VuTLPs were induced in the first hours in response to the MI_CABMV or MI_CPSMV. There are previous reports in the literature about the interaction between viral proteins and TLPs. Kim et al. (2005) observed that a Nicotiana tabacum (NtTLP1) TLP specifically interacted with CMV (Cucumber mosaic virus) - ‘movement-related protein’, ‘movement protein’ and ‘coat protein’ - during a yeast two-hybrid screen experiment. However, the determination of the possible role of TLP-virus interaction requires deepened studies.

When comparing the transcriptional orchestration in response to MI_viral, the results indicate a functional specialization of VuTLPs in cowpea submitted to MI_CABMV and MI_CPSMV, since only specific representatives are induced, in detriment of others, in response to inoculation by a given virus / injury.

The results for radicular dehydration assay indicate that some TLPs are recruited in the first minutes after stress imposition, even though a higher number of transcripts were detected at 150 min. Previous studies did not address the global transcriptional orchestration of TLP proteins and referred to individual case studies. An example is given by Jung et al. (2005) that reported a new TLP isoform isolated from carrot, regarding a drought-specific, ABA-independent, non-organ-specific, and non-developmental-stage-specific protein.

The comparison between the orchestration in response to both viruses (in leaves) and radicular dehydration (roots) revealed potential VuTLPs expressed in both assays, indicating a crosstalk role for TLPs under different stress types and tissues. Besides, there is no overall conservation in VuTLPs orchestration towards different viral isolates and radicular dehydration, indicating specific TLP associated with the response of each virus type and radicular dehydration.

Signalization response in plants submitted to abiotic and biotic stressors can induce separate and overlapping sets of genes, leading to the expression of distinct as well as common components (Zhu et al. 2014; Mellacheruvu et al. 2016). These separate pathways show nodal points where they converge and crosstalk to optimize the various defense responses (Xiao et al. 2013), resulting in shared stress mitigation strategy by combined morpho-physiological processes and molecular responses (Pandey, Ramegowda, Senthil-Kumar, 2015). The identification of crosstalk between biotic and abiotic signaling pathways has been crucial for envisaging and strengthening our understanding of the regulation of plants response against combined stresses. Genes such as those coding for TLPs that impart abiotic as well as biotic stress tolerance are highlighted in some studies in order to understand the mechanism underlying tolerance to dual stresses.

Singh et al. (2013) analyzed an Arachis diogoi (AdTLP) TLP induced under fungal infection (Phaeoisariopsis personata). The cDNA encoding AdTLP was cloned using RACE-PCR and used to transform tobacco plants. Overexpression of AdTLP resulted in increased resistance to pathogenic fungi and tolerance to abiotic stresses (high salinity and oxidative stress). In addition, transgenic plants also exhibited a higher level of transcription of genes PR1a, PI-I and PI-II compared to WT. Such genes are associated with plant pathogen defense mechanisms. Chowdhury, Basu and Kundu (2017) showed that transgenic lines of sesame overexpressing an osmotin-like protein (SindOLP) presented tolerance against abiotic stresses (drought and high salinity) and resistance against the fungus Macrophomina phaseolina. Overexpression of SindOLP resulted in the induction of three genes [superoxide dismutase (SiSOD), cysteine protease inhibitor (SiCysPI), and glutathione-S-transferase (SiGST)] which encode enzymes for the elimination of reactive oxygen species (ROS), indicating that SindOLP participates in the ROS regulation, which is common to both stresses addressed.

Considering the results of qPCR, it is evident that TLPs are recruited in the first minutes after MI_CABMV or MI_CPSMV treatment. Even though both treatments had the injury in common, different TLPs were recruited after each virus inoculation, indicating some specificity that deserves additional evaluations.

Some recently published papers have shown the induction of TLPs after virus inoculation. Madroñero et al. (2018) analyzed changes in the papaya transcriptome in response to the inoculation with PMeV complex (Papaya Meleira Virus, PMeV, and Papaya Meleira Virus 2, PMeV2). To evaluate the effects of SA signaling on PMeV complex load, Carica papaya seedlings were exposed to this hormone and inoculated with both viruses. The SA-treated plants showed increased PR1, PR5 and CHIA levels confirming SA signaling activation. The same plants presented reduced PMeV complex load, suggesting that SA signaling plays a role in C. papaya tolerance to PSD before flowering.

Kappagantu et al. (2017) studied the effect of the host response of hop plants (Humulus lupulus) to Hop stunt viroid (HSVd) on the development of hop powdery mildew (Podospheara macularis; HPM). Transcriptome analysis followed by qPCR analysis showed that transcript levels of PR-genes such as PR-1, chitinase, and TLPs were induced in leaves infected with HPM alone. The response of these genes to HPM was significantly down-regulated in leaves with HSVd-HPM mixed infection. These results confirm that HSVd alters host metabolism, physiology, and plant defense responses.

The results of qPCR show that – despite being traditionally responsive to biotic stresses – TLPs have been reported as recruited under abiotic stresses (Zhao et al. 2024; Muoki et al. 2021; Jung et al. 2005; Ruperti et al. 2002; Petre et al. 2011). In the present work, a higher number of VuTLPs were induced at 150 min after root dehydration. Zhao et al. (2024) identified and characterized 10 TLP genes from the faba bean transcriptome. Two genes (VfTLP4-3 and VfTLP5) displayed high expression levels in the drought-tolerant cultivar when exposed to drought stress. In turn, Misra et al. (2016) identified and characterized a new basil (Ocimum basilicum) TLP (ObTLP1) from ESTs recovered after MeJA treatment. The ectopic expression of ObTLP1 in A. thaliana provided fungal resistance (Scleretonia sclerotiorum and B. cinerea) and tolerance to abiotic stresses (drought and high salinity).

Despite the evidence of the present work and previous studies, the mechanism of action of such proteins against abiotic stress remain largely unknown, although investigations have confirmed the roles of several TLP genes in cold temperature, salt and drought stresses (Jung et al. 2005; Muoki et al. 2012; Singh et al. 2013; Misra et al. 2016; Zhao et al. 2024). Zhao et al. (2024) showed ectopically expression of two genes from the faba bean (VfTLP4-3 and VfTLP5) in tobacco increased drought tolerance. This enhancement in drought tolerance was accompanied by changes in the activity of peroxidase (POD) and superoxide dismutase (SOD), as well as the levels of malondialdehyde (MDA) and proline content. VfTLP4-3 appeared to primarily regulate chitin metabolism, while VfTLP5 was associated with the regulation of proline synthesis. The work of Chowdhury, Basu and Kundu (2017) with sesame, indicated a role in the induction of three genes [superoxide dismutase (SiSOD), cysteine protease inhibitor (SiCysPI), glutathione-S-transferase (SiGST)] which encode enzymes for the elimination of reactive oxygen species (ROS). These results suggest that SindOLP participates in the regulation of the ROS pathway, which is common to both types of stresses.

5. Conclusion

The present study describes the structural diversity and gene expression profile of thaumatin-like proteins (VuTLPs) in cowpea under biotic and abiotic stress conditions. Identification of 34 loci encoding VuTLPs along with the detection of both segmental and tandem duplication events underlines the complexity of this gene family in cowpea and its evolutionary dynamics. The prevalence of purifying selection across VuTLPs suggests that these proteins have conserved functions important for plant adaptation to stress. The obtained results showed that different stressors provoke distinct expression profiles in VuTLPs and that the most rapid induction generally occurs within minutes following the onset of stress. Delayed expression in response to root dehydration underlines the possibility of functional specialization associated with certain types of unfavorable conditions. This is further supported by the analysis of promoters linking VuTLPs with various transcription factors, including bHLH, Dof-type, and MYB-related factors, indicating multifaceted participation in response mechanisms. Collectively, these analyses provide an overview of VuTLPs and suggest them as potential targets for genetic engineering strategies in cowpea. Towards this, future research studies should henceforth be directed to the elucidation of VuTLP-mediated molecular pathways and interaction with other stress-related proteins for fullest biotechnological applications in crop improvement.