The complete genomic sequence of the pea plant was sourced from the National Center for Biotechnology Information (NCBI) database. A hidden Markov model (HMM) was crafted using a variety of kinesin (KIN) protein sequences from different species, which are accessible on the NCBI (GCA_024323335.1) platform. The HMMER software package (version 3.3.2) was then applied to perform searches against the local protein database of the pea plant, with an E-value threshold set at 10^-5 to filter results (Supplementary Table S1). This process resulted in the identification of a preliminary list of KIN candidate genes, which was subsequently refined by eliminating any duplicates. This method has been reported in the literature (Potter et al., 2018). In rice, Arabidopsis, and soybean, the HMMER model construction was also used (Guo et al., 2008; Chen et al., 2014).

To substantiate the authenticity of these candidate sequences, the NCBI Conserved Domain Database (CDD) and the SMART tool were employed. These tools were used with an E-value threshold of 10^-5 to filter the sequences based on the presence of the characteristic KIN domain. The KIN genes that passed this validation were systematically renamed according to their respective chromosomal locations within the pea genome.

For the prediction of subcellular localization, the Cell-PLoc software (version 2.0) was engaged. Additionally, the ExPASy ProtParam tool was leveraged to forecast various physicochemical attributes of the proteins, encompassing molecular weight (MW), isoelectric point (pI), instability index, and grand average of hydropathicity (GRAVY).

To elucidate the phylogenetic connections among the kinesin (KIN) genes of the pea plant, we conducted multiple sequence alignments involving the identified pea KIN proteins and those from Manihot esculenta, Populus trichocarpa, and A. thaliana. The sequence alignment was executed using the ClustalW tool, and the resultant data were utilized to construct a maximum likelihood (ML) phylogenetic tree with IQ-TREE 2 software, applying a bootstrap value of 1,000 to ensure the robustness of the tree topology.

The chromosomal positions and the exon-intron structures of the KIN genes were extracted from the genome annotation files. These files were sourced from the Ensembl Plants database (http://plants.ensembl.org/index.html). The physical locations of the KIN genes on the chromosomes were visualized using the MapGene2Chromosome v2.0 (http://mg2c.iask.in/mg2c_v2.0/). The exon-intron structures were graphically represented with the aid of the Gene Structure Display Server (GSDS2.0) (http://gsds.cbi.pku.edu.cn/).

Furthermore, to identify the conserved domains within each KIN protein sequence, we employed the MEME tool (http://meme-suite.org/). The search for motifs was configured to identify up to 10 distinct motifs, with all other parameters left at their default settings to ensure a thorough and unbiased search for conserved elements.

The genomic coordinates and structural attributes of genes on the chromosomes of pea were extracted from the GFF3 file, which was sourced from the pea genome database (htKIN://www.peagdb.com/download/). Subsequently, the MapGene2Chrom web v2 (http://mg2c.iask.in/mg2c_v2.0/) was applied to graphically represent the PsKIN genes on their respective chromosomes. The analysis of synteny among the KIN genes across pea, A. thaliana, and Glycine max was conducted using the Multiple Collinearity Scan Toolkit (MCScanX) and TBtools. The genomic datasets for A. thaliana and G. max were obtained from the Phytozome12 database (htKIN://phytozome.jgi.doe.gov/pz/portal.html).

The analytical capabilities of MCScanX were leveraged to determine the Ka (non-synonymous substitution rate) and Ks (synonymous substitution rate) for gene pairs involved in both segmental and tandem duplications. The Ks values were used to calculate the dates of duplication events (T) according to the following equation: T = Ks/2λ, with λ = 1.5 × 10⁻⁸ s for dicots (Bian et al., 2020). The Ka/Ks value was further employed to identify the selection mode of the PsKINs. STRING (htKIN://string-db.org/) was used to construct the functional interaction network of the proteins.

The propagation of the ‘Zhewan No.1’ (ZW1) TM-1 seeds of variety was initiated in a substrate composed of nutrient-enriched soil within an artificially regulated greenhouse. The conditions within the greenhouse were tailored to provide a 16 h period of light exposure at a steady temperature of 27°C, followed by an 8 h interval of darkness maintained at 22°C, with the atmospheric humidity consistently maintained between 60% and 80%. Peas seedlings of the TM-1 variety, which had reached a developmental stage characterized by the presence of three leaves and demonstrated signs of robust health and stable growth, were subsequently introduced to treatments designed to induce drought and salt stress. For the induction of salt stress, the root systems of the seedlings were treated with an irrigation solution containing a concentration of 300 mM sodium chloride (NaCl). In contrast, to mimic the effects of drought stress, a 20% solution of polyethylene glycol 6,000 (PEG6000) was applied to the root systems. Following the application of these stressors, leaf samples were harvested at several time points post-treatment: immediately post-treatment (designated as 0 h), as well as at 3 h and 24 h intervals. Upon collection, each leaf sample was rapidly immersed in liquid nitrogen to achieve rapid freezing, and then securely stored at a temperature of −80°C for the purpose of RNA extraction in subsequent experimental procedures.

Post-treatment at 3 h and 24 h time points, the fresh pea seedling leaves from the untreated control group and those subjected to 300 mmol NaCl for both time intervals, as well as the leaves treated with PEG6000, were flash-frozen in liquid nitrogen and dispatched to Tianjin Jizhi Gene Technology Co., Ltd. For further analysis. In order to validate the uniformity across our experiments, triplicate biological samples were analyzed. RNA integrity was determined by agarose gel electrophoresis, revealing clear 28S and 18 S rRNA bands in the total RNA, with the 28 S intensity of band being roughly double that of the 18 S band. For the preparation of the RNA-seq library, a multi-step protocol was adhered to, encompassing total RNA isolation, mRNA purification, fragmentation, cDNA generation, end-repair of double-stranded cDNA, adapter attachment, and PCR amplification. Additionally, the RNA-seq dataset was subjected to Weighted Gene Co-expression Network Analysis (WGCNA) to assess and secure the precision and dependability of the transcriptome data.

The resulting raw sequencing reads underwent a stringent quality control process and were subsequently cleaned. These high-quality reads were aligned to the reference genome, accessible via the provided link, to facilitate the assembly of transcripts and the quantification of gene expression levels. The quality of the RNA-seq alignments was meticulously evaluated. The subsequent analysis of the transcriptome was performed utilizing a proprietary cloud-based platform developed by Tianjin Jizhi Gene Technology Co., Ltd., ensuring a comprehensive and in-depth examination of the data.

In our research, we initiated the investigation by employing a hidden Markov model (HMM) specific to plant kinesin (KIN) proteins, utilizing this model to conduct a comprehensive search across the entire proteome of pea. Subsequently, we proceeded to authenticate the credibility of the identified candidate KIN proteins. A total of 105 KIN proteins, characterized by the presence of a conserved KIN domain, were identified and designated as constituents of the pea KIN family, with each member named PsKIN1 through PsKIN105, corresponding to their chromosomal positions (Figure 1; Supplementary Table S2).

Further analysis was conducted to elucidate the fundamental characteristics of the PsKIN genes, encompassing parameters such as the length of the protein in amino acids (aa), molecular weight (MW), isoelectric point (PI), instability index, and grand average of hydropathy (GRAVY). The length of the PsKIN proteins varied significantly, extending from 70 aa for PsKIN60 to a substantial 2,980 aa for PsKIN101, with an overall mean length of 541.2 aa. The molecular weight of these proteins also exhibited a broad range, from 7.66 kD for PsKIN60 to a notably higher value of 341.67 kD for PsKIN101, with an average MW of 76.50 kD. Notably, both the length and MW of PsKIN101 exceeded those of the other PsKIN genes by a considerable margin. Additionally, the isoelectric points of the PsKIN genes spanned from 4.5 to 9.94, the instability indices ranged from 34.84 to 58.89, and the GRAVY values varied from −0.867 to 0.088.

To elucidate the evolutionary connections among the KIN gene family members in pea (comprising 105 PsKINs), Arabidopsis (encompassing 65 AtKINs), and G. max (consisting of 150 GmKINs), we developed a maximum likelihood (ML) phylogenetic framework. This framework was utilized to delineate the phylogenetic affinities among the species in question (Figure 2). Within the ML phylogenetic tree, the 320 KIN genes were categorized into seven distinct subfamilies, following the taxonomic and topological criteria established by prior research. The distribution of the 105 PsKINs across these subfamilies was found to be non-uniform, as detailed in Supplementary (Supplementary Table S2). Among the subfamilies, Subfamily III was identified as the most populous, housing 25 PsKIN genes. In contrast, Subfamily I was the least represented, with only 5 PsKIN genes. This differential distribution underscores the diversity and evolutionary complexity within the KIN gene family across the examined species.

In our investigation, we conducted a comprehensive analysis of the gene structure and the conserved motifs present within the 105 PsKIN proteins, correlating these findings with their phylogenetic relationships (Figures 3A–D). Utilizing the MEME online tool, we identified 10 conserved motifs (motifs 1–10; Figure 3B) across the PsKIN proteins. The analysis highlighted variations in the composition of these motifs, with a general trend showing that the motif distribution aligns well with the classification into subfamilies. Notably, motifs 5, 6, and 7 emerged as the predominant motifs across all PsKIN genes.

Concurrently, we examined the gene domains of the PsKINs, presenting a visual representation of the conserved KIN domains mapped directly onto the genetic structure of the pea (Figure 3C). Our findings indicated the presence of 12 distinct types of KIN domains, with the majority of genes anticipated to encompass the “Kinesin motor domain”.

Further analysis of the exon count within the 105 PsKINs, as per the pea GFF annotation file, revealed a range from 1 to 36 exons (Figure 3D). A significant number of PsKINs were found to possess more than 10 coding sequences (CDSs). PsKINs from related subfamilies displayed a tendency towards similar exon/intron structures, indicating a strong correlation between phylogenetic relationships and gene structure within the gene family. This correlation underscores the evolutionary significance of gene structure in shaping the diversity and function of the KIN gene family in pea.

Cis-regulatory elements (CREs), which are non-coding DNA sequences located within the promoter region of genes, play a pivotal role in gene expression and are integral to the regulation of a multitude of biological processes. In our study, we harnessed the 2000 base pair (bp) promoter sequences upstream of the identified PsKIN genes to forecast the presence of CREs, utilizing the PlantCARE database as a predictive tool. A total of 2,111 CREs sites, deemed representative, were culled from the predictive outcomes and are presented in Figure 4A. Within this set, elements associated with growth and developmental processes were found to be the most prevalent, closely trailed by those linked to hormone responsiveness, as illustrated in Figure 4B.

This analysis underscores the complexity of gene regulation within the PsKIN family, highlighting the significance of promoter-associated CREs in modulating gene expression in response to various developmental cues and hormonal signals. The abundance of growth and development elements, as well as hormone responsiveness elements, suggests a dynamic interplay between these regulatory factors and the expression of PsKIN genes, potentially influencing their roles in plant physiology and adaptation.

The CREs on the putative promoter of the PsKINs. (a) Distribution of CREs identified in the 2000 bp upstream promoter region of PsKINs; (b) The number of CREs on putative promoters of PsKINs. The numbers in the heatmap represent the quantity of elements.

To delve into the genomic expansion mechanisms of the PsKINs gene family within pea, an examination of the syntenic relationships among the PsKINs was conducted. Out of the 105 PsKINs, 37 exhibited syntenic relationships, with 21 of these forming syntenic pairs that have undergone segmental duplication events within the pea genome, as detailed in Supplementary Table S4. The ratio of the rate of non-synonymous substitutions to the rate of synonymous substitutions, denoted as Ka/Ks, serves as an indicator of selective pressures acting on genes post-duplication. For the 21 segmental duplication gene pairs, the Ka/Ks ratios were found to be less than 1, suggesting that these pairs have experienced strong negative selection throughout their evolutionary trajectory.

In an effort to enhance our comprehension of genetic divergence, gene duplication, and evolutionary patterns between the KIN gene families of P. sativum, A. thaliana, and G. max, syntenic relationships were scrutinized to pinpoint orthologous KIN genes among these species using the MCScanX tool. A total of 56 and 184 pairs of orthologous KIN genes were identified across three comparative analyses (P. sativum versus A. thaliana, P. sativum versus G. max), as depicted in Figure 5. It is noteworthy that certain genes exhibit multiple homologous relationships across the different species. Notably, a significant concentration of syntenic genes was observed on chromosomes Ps5, Ps6, and Ps7.

To gain a deeper insight into the biological roles and the intricate regulatory networks associated with PsKINs, the potential protein-protein interactions (PPIs) among these proteins were forecasted employing an orthology-based approach. The outcomes revealed that 43 of the PsKIN proteins shared orthologous connections with their counterparts in Arabidopsis. Intriguingly, this analysis identified the presence of both physical interactions and co-expression patterns among these 43 PsKIN proteins (Figure 6, Supplementary Table S5). It is hypothesized that these interactions are instrumental in orchestrating plant growth and development through the modulation of protein networks.

To gain a comprehensive understanding of the expression profiles of the PSKIN gene family across various tissues and to anticipate their possible biological roles, an expression heatmap was developed, encompassing 38 PSKIN genes (Figure 7). These were meticulously picked from a collection of 105 genes, showing distinct patterns in a range of tissues including primary and lateral roots, nodules, young and mature stems, tendrils, young and mature leaves, sepals, petals, pods, and seeds.

We identified that PSKIN76 is markedly expressed in the roots, with PSKIN96 showing elevated levels specifically in lateral roots. Conversely, expression peaks of PSKIN82 in young stems, petals, and pods. Additionally, there is a notable upregulation of PSKIN103 in young stems and seeds, while PSKIN17 and PSKIN98 are highly expressed in seeds.

Our in-depth analysis indicated that 20 genes are predominantly expressed in seeds, hinting at their involvement in seed maturation and dormancy. Meanwhile, 11 genes displayed high expression specificity in pod tissues, potentially linked to pod development and seed safeguarding mechanisms.

To explore the reaction of the KIN gene family under drought conditions, pea seedlings were placed in an environment with a 20% concentration PEG6000 solution. Leaf tissues were harvested at 3 h and 24 h intervals post-treatment for subsequent RNA extraction and transcriptome profiling. As illustrated in eight of picture, a notable fluctuation in the expression levels of 47 PsKIN genes was discerned, including PsKIN8, PsKIN11, PsKIN36, PsKIN37, PsKIN54, PsKIN102, and PsKIN104, following exposure to PEG6000 (Figure 8). This biphasic response, characterized by an initial upregulation followed by a downregulation, could be indicative of acute reaction of the plant to water scarcity and the subsequent acclimation phase. Preliminary functional analyses hint that these genes might partake in signal transduction and metabolic pathways that are instrumental in the plant’s drought stress management.

Conversely, a subset of genes, namely, PsKIN26, PsKIN47, PsKIN54, and PsKIN64, exhibited a sustained decrease in expression, potentially highlighting their involvement in long-term stress endurance mechanisms. This sustained response may be crucial for the preservation of cellular hydration and the reinforcement of antioxidant systems.

Expanding our inquiry to the realm of saline stress, seedlings were treated with escalating concentrations of NaCl, 100, 200, and 300 mM. Post-treatment analysis revealed pronounced alterations in the expression profiles of numerous PsKIN genes, with PsKIN47, PsKIN51, and PsKIN90 standing out as the most significantly affected. These changes are presumed to be intimately linked to the mechanisms of the plant of acclimatization to high salinity, possibly implicating these genes in osmotic adjustment and ionic equilibrium.