The sequence IDs of Arabidopsis and rice RPK1 genes were acquired from the available literature (Shi et al., 2014; Zou et al., 2014). These sequences were retrieved from Ensembl plants and NCBI, which were then used as queries for the Basic-Local Alignment Search tool (BlastP and BlastN) against IWGSC (INSDC Assembly GCA_900519105.1 July 2018 database version 106.4), NCBI (https://www.ncbi.nlm.nih.gov/) and Ensembl plants (plantshttps://plants.ensembl.org/index.html) for T. aestivum. For all of the candidate RPK1 genes, the kinase domain presence was substantiated with Pfam (http://pfam.sanger.ac.uk), and by SMART (http://smart.embl-heidelberg.de/) (Letunic and Bork, 2018) databases. The sequences in which the kinase domain was absent were removed (Supplementary Table S1). In silico based putative protein information of RPK1 genes (physio-chemical) was analyzed through the Protparam (https://web.expasy.org/protparam/) tool. The subcellular localization of RPK1 proteins was predicted via Plant-mSubP and pLoc-mPlant (http://bioinfo.usu.edu/Plant-mSubP/ ; http://www.jci-bioinfo.cn/pLoc-mPlant/ ) (Cheng et al., 2017; Sahu et al., 2020).

The chromosomal locations of all candidate RPK1 genes in T. aestivum were acquired from Ensembl (http://plants.ensembl.org/Triticum_aestivum/Info/Index). The gene map of TaRPK1 genes was drawn with the help of MapChart and confirmed through TBtools.

To retrieve the RPK1 protein sequences, the amino acid sequences of 15 TaRPK1 members were used as queries to blast (BLASTP) against the Triticum turgidium, Triticum dicoccoides, Titicum urartu, Triticum speltoides, Aegilops tauschii, Hordeum vulgare, Arabidopsis thaliana, and different species of Oryza (rufipigon, japonica, indica, and glaberrima). The sequences with > 60% identities were retrieved from Ensembl (http://plants.ensembl.org). The phylogenetic trees were made by means of MEGA-X software with NJ (neighbor-joining method) (Kumar et al., 1994). The parameter Poisson model and pairwise deletion were used with replicates of 1,000 bootstraps for assessment of node significance.

The number of exons and introns was predicted by the gene structure display server (GSDS, http://gsds.cbi.pku.edu.cn/) and the genomic sequences and coding sequences were aligned using ClustalW. Conserved motifs in RPK1 proteins of T. aestivum were analyzed using MEME, a multiple-EM for motif elicitation program (http://meme-suite.org/tools/meme) (Bailey et al., 2009). The execution of MEME search was done with default parameters apart from motif maximum number, which was set to 10, and optimum motif width of ≥6 and ≤200 was selected.

The analysis of TaRPK1 gene ontology was performed by TaRPK1 protein sequences via the online gProfiler tool (https://biit.cs.ut.ee/gprofiler/gost) with default parameters (Raudvere et al., 2019).

The miRNA prediction was performed as mentioned formerly (Yan et al., 2019). The TaRPK1 sequences were submitted for potential miRNA prediction through a search against the available wheat miRNA reference by means of the psRNATarget Server (https://www.zhaolab.org/psRNATarget/), using default settings (Dai and Zhao, 2011). The visualization of the interaction network of the predicted miRNA with their corresponding TaRPK1 target genes was done by Cytoscape software (https://cytoscape.org/) with default settings (Shannon et al., 2003).

The sequence size of 2 kb in the upstream region were dug out from all TaRPK1 genes of T. aestivum that acted as promoters for the regulatory cis acting elements prediction through the PlantCare (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) database (Lescot et al., 2002).

The GFF3 files and proteomes of Triticum aestivum and its ancestors, including Aegilopes tauschii, Triticum spelta, Triticum turgidum, and Triticum dicoccoides, were used from the Ensembl Plants database for collinearity prediction via the MCScanX algorithm (Wang et al., 2012). Synteny scrutiny of RPK1 family members was performed via Tbtools (Chen et al., 2018).

The TaRPK1 protein structures were modeled via amino acid sequence using the SWISS-MODEL database (https://www.swissmodel.expasy.org) (Biasini et al., 2014), and for visualization of 3D structure Pymol software (https://pymol.org/2/) was applied. The verification and validation of the predicted 3D structures of TaRPK1 proteins were assessed using the Ramachandran Plot—Zlab, (https://zlab.umassmed.edu/bu/rama/) (Anderson et al., 2005).

In silico expression analysis was performed using the wheat-expression browser (www.wheat-expression.com) at different wheat stages (Kaur et al., 2017). The data were unruffled in the course of developing seedling, vegetative, and reproductive stages from different organs of wheat such as roots, leaf sheath, leaf blade, shoot, spike, and grain. The heatmap was then created from the composed data, based on the expression values of genes (in TPM) by means of Tbtool.

For interaction network studies, String (https://string-db.org/) was used by selecting Triticum aestivum as a platform species. For visualization of the molecular library, Cytoscape was used. Correlation coefficients on the basis of verities, treatments, and tissues were calculated in R 3.4.0. These coefficients indicate the degree of association among the terms and provide linkages among the TaRPK1 members.

Previously, Pakistan-13, Galaxy, and Shafaq were studied under drought stress and categorized as drought tolerant and susceptible varieties, respectively (Shabbir et al., 2015; Ulfat et al., 2017; Ahmad et al., 2021; Wasaya et al., 2021; Iqra et al., 2022). So, seeds of these varieties were obtained and sown under controlled glass-house conditions at the National Institute for Genomics and Advanced Biotechnology (NIGAB), National Agriculture Research Center (NARC), Islamabad, Pakistan. After 2 weeks of sowing (seedling stage), the roots and leaves tissues were collected. At growth stage 8 (tillering stage), roots, stems, and leaf tissues were collected. At the grain filling stage (14 days after flowering), sampling for roots, stems, leaf, and grains was done (Hyles et al., 2020). For expression profile analysis under drought stress, seeds of selected varieties were first surface sterilized with sodium hypochlorite followed by three washings, then soaked in distilled water in a growth chamber (16 h light/8 h dark cycle at 22°C). After 2 weeks, young seedlings were treated with 20% polyethylene glycol (PEG) 6,000 (v/v). The root and leaf tissues of seedlings were harvested after 12 h of exposure to stress conditions. All the samples were collected in three replicates, and samples were frozen immediately in liquid nitrogen, and placed in −80°C storage for RNA extraction.

Approximately 100 mg of tissues were taken for total RNA extraction using an RNA mini kit (Cat # 12183018A, Invitrogen, Thermo Fischer Scientific) followed by the manufacturer’s instructions. Through agarose gel electrophoresis, the quality and concentration of RNA were determined, followed by optical density measurement through a spectrophotometer. With the help of the RT Prime-Script Reagent Kit, the cDNA was made from 1 ug of RNA. Specific primers were designed for TaRPK1 genes manually, followed by confirmation via NCBI Primer Blast software (http://www.ncbi.nlm.gov/tools/primer-blast), provided in Supplemental Table S2. The qRT-PCR was accomplished with SYBR Green I (Roche) Master Mix. Wheat β-Actin was used as a control reference gene. Three independent biological replicates were analyzed for each sample. The values were means and standard deviations (SD) were calculated from biological replicates. The relative expression levels of each gene were studied by means of 2^−∆∆Ct (Schmittgen and Livak, 2008).

A set of 15 candidate RPK1 genes were retrieved from Triticum aestivum based on BlastP and BlastN. A domain search by the SMART tool with the corresponding RPK1 candidate amino acid sequences confirmed the S_TKc Domain (SM00220). Thus, a total of 15 TaRPK1 with complete structures were analyzed in T. aestivum (Table 1). Subsequent sequence identification of 15 TaRPK1 showed the protein length of 609–1,123 amino acids and a molecular mass ranged from 67–120 kDa. The iso-electric points (PI) of these proteins were 6–9. The Instability Index (II) ranged from 28.37–52.17, the Aliphatic Index (AI) was 87.93–106.38, and the grand-average of hydropathicity (GRAVY) −0.032–0.138. The instability index of group I was less than 40, representing stable proteins, whereas proteins of groups II and III showed instability index values of more than 40, indicating unstable proteins. The AI signified that all of the TaRPK1 proteins are thermally stable. The GRAVY indicated TaRPK1 proteins to be hydrophilic proteins except for TaRPK10, TaRPK11, and TaRPK12, which showed a value less than zero, representing them as hydrophobic proteins. The sub-cellular localizations of the TaRPK1 were anticipated, which showed that all the TaRPK1 were localized to the cell membrane (Table 2).

A detailed protein alignment of structural predictions showed that all TaRPK1 proteins are composed of the leucine-rich repeat N terminal (LRRNT_2) domain, leucine-rich repeat (LRR) domains, transmembrane domain I, and a serine–threonine kinase (S_TKc) domain. However, the LRR domains were missing in TaRPK1, TaRPK2, and TaRPK3 sequences (Supplemental Figure S1).

The physical location of RPK genes in T. aestivum, to the corresponding chromosomes, is shown in Figure 1. A total of 15 RPK genes were mapped on 9 out of 21 chromosomes in wheat. The genes were mainly mapped on chromosomes 1, 2, and 3 on the respective A, B, and D genomes. No RPK genes were found on the rest of the chromosomes.

Of the 15 identified TaRPK1 genes in this study in Triticum aestivum, two RPK genes from Arabidopsis thaliana, 16 RPK genes from rice, seven RPK genes from Triticum dicoccoides, three RPK genes from Triticum urata, seven RPK genes from Triticum turgidium, four RPK genes from Aegilops tauschii, 11 RPK genes from Triticum speltoides, and four RPK genes from Hordeum vulgare were used to construct a neighbor-joining based tree with MEGA X software in order to study the evolutionary relationships (Figure 2). The phylogenetic tree generated on the basis of similarities with protein sequences distributed RPK members into four main groups, with TaRPK1 members in three groups. Overall group I possessed nine TaRPK1 members (TaRPK1-9), that were closely associated with RPK members of rice. Group II (TaRPK10-12) and Group III (TaRPK13-15) exhibited three TaRPK1 members each, that exhibited close association with Triticum turgidium, Triticum speltoides, Aegilops tauschii and Triticum dicoccoides.

The intron–exon number and arrangements of the RPK1 members were envisaged through comparing the coding sequence with the genomic DNA sequence. All of the TaRPK1 genes in group I consisted of 16–18 introns, except for the groups II and III that contained 0 and 1 intron (Figure 3A). Furthermore, the conserved motifs within TaRPK1 proteins were predicted by online MEME software. Ten conserved motifs (1-10) were analyzed (Figure 3B and Supplementary Figure S2). The motifs 1, 3, 4, 7, 8, and 10 were present in all of the RPK1 sequences. However, group II did not display motifs 3 and 10, and the motif three was also missing in group III sequences.

GO annotation analysis was conducted for the functional analysis of RPK1 genes. In-silico functional prediction was performed, and the results displayed two types of processes involved, that is, molecular processes (MPs) and biological processes (BPs) (Figure 4 and Supplementary Table S3). Biological processes indicate that RPK members are involved actively in various metabolic processes. The molecular processes suggested the RPK1 member’s catalytic activity. Such outcomes clearly denote RPK1 genes’ significant role in growth and development via modulation of molecular and biological processes.

We also identified putative 18 miRNAs targeting TaRPK1 genes for the generation of interaction networks by Cytoscape software in order to better understand the underlying miRNA mechanism involved in the modulation of TaRPK1 genes (Figure 5 and Supplementary Table S4). In the connection distribution and regulation network, TaRPK1, TaRPK2, and TaRPK3 were found targeted by single miRNAs, which are tae-miR9782, tae-miR9776, and tae-miR1122c-3p, respectively. TaRPK10 and TaRPK11 are the most targeted RPK1 wheat genes by tae-miR1134, tae-miR9774, tae-miR9661-5p, tae-miR9664-3p and tae-miR9777 targeting TaRPK10, and tae-miR9774, tae-miR9777, tae-miR9664-3p, tae-miR395a and tae-miR9661-5p targeting TaRPK11 genes. However, no miRNA was found targeting TaRPK13, TaRPK14, and TaRPK15 genes.

The promoter regions contain cis-modulatory elements which are critical for the binding of transcription factors for transcription initiation, which has an essential function in the expression of genes. The promoter regions of RPK1 members were used for the cis-regulatory element prediction (Figure 6A). The results indicated that the cis-regulatory elements can be distributed into several categories, such as hormone related elements, light-related elements, developmental responsive elements, abiotic stress responsive elements, promoter-related motifs, and other motifs. Amid them, the elements chiefly present were associated with photoreaction, hormone responsiveness, and abiotic stress-related motifs. The photoreaction responsive cis-regulatory elements included ACE, AE-Box, ATCT, G-Box, GATA, GT1, SP1, AT1, Box 4, Box II, I-Box, TCT, GA, L-Box, TCCC, and ATC motif. The most abundant light-responsive elements were found in TaRPK11 and TaRPK13, which had 17 and 12 members, respectively. Hormone responsive elements were also copiously present in the RPK1 promoter, mostly comprising abscisic acid response elements. The three extensively distributed cis elements were related to abiotic stress response, among which drought responsive elements were profuse. Other elements correlated to abiotic stress were also identified.

In order to understand the evolutionary relationship and origin of Triticum aestivum (tr) with Triticum turgidum (tg), Aegilops tauschii (at), Triticum speltoides (ts) and Tritium dicoccoides (td), a comparative synteny scrutiny of RPK protein sequences was performed. The proteins were closely related among five species and exhibited significant similarity in analysis of evolutionary correlation. It was observed that the TaRPK1 genes of T. aestivum have similar origins of evolution to other Triticum species (Figure 6B and Supplementary Table S5).

Three-dimensional (3D) structures of TaRPK1 proteins were predicted by using SWISS_MODEL online computational software. 3D structures of target proteins were anticipated based on homology modeling. The SWISS MODEL predicted 15 successful models of TaRPK1 proteins with at least 30% identity to the template (4mn8.1. A, 5hyx.1. A, 5xkj.1. C, 6mOu.1. A, 4mna.1. A, 4oh4.1. A, 6cth.1. A, 7brc.1. A, and 5tos.1. A) that was a widely recognized threshold for effective modeling (Xiang, 2006). However, TaRPK2 and TaRPK3 showed sequence identity of 27.84% and 29.47%, respectively, with the template, which was less than 30%. The highest sequence identity of 45% with the template was observed by TaRPK4, TaRPK5 and TaRPK6 (Figure 7). The verification and validation of the predicted 3D structure of TaRPK1 were assessed via Ramachandran Plots (Anderson et al., 2005) that validated the backbone diahedral angles of the targeted protein. The Ramachandran plot assessment showed that 92–98% of the regions of TaRPK1 protein showed highly favorable regions, which indicates the stability and good quality of the predicted protein structure (Supplemental Table S6).

The data of RNA-seq for all of the 15 RPK sequences were obtained from online database. A heatmap was generated showing expression levels of RPK members at different stages, namely seedling stage, vegetative stage, and reproductive stage (Figure 8 and Supplementary Table S7) and in various organs (root, leaf, shoot, spike, and grain) of wheat. The highest expression of TaRPK1 members was observed in root tissues compared to other tissues. TaRPK1, TaRPK2, and TaRPK3 exhibited the highest expression patterns in roots at seedling, vegetative, and reproductive stages. Higher to moderate expression was observed in grain at the developing reproductive stage by TaRPK13 and TaRPK14, respectively. Spikes, leaves, and shoots showed moderate to low expression in all of the TaRPK1 members in wheat.

The TaRPK1 gene expression was determined in drought-tolerant (Pakistan 13 and Galaxy) and drought-susceptible (Shafaq) varieties under normal growth conditions in order to get a baseline expression profile. The expression pattern in all of the three varieties was examined in various developmental stages, including seedling stage, tillering stage, and heading stage and in different tissues such as root, stem, leaf, and grain (Figure 9). The TaRPK1, TaRPK2, and TaRPK3 showed significant expression in the roots at the heading and seedling stages of the Pakistan-13 and Galaxy varieties. The TaRPK13 exhibited higher expression in grain tissues of all varieties compared to other TaRPK1. The TaRPK1 genes displayed higher expression in roots whereas they showed less expression in leaves and stems compared to the grain and root expression in developmental stages. Our results indicated that TaRPK1 genes had similar expression patterns in both Pakistan 13 and Galaxy varieties, unlike the Shafaq variety. The higher expression of TaRPK1 genes was observed in the heading > seedling > tillering stages in Pakistan 13, Galaxy, and Shafaq varieties. Overall, TaRPK1 exhibited significant expression in root tissues compared to leaf, shoot, and grain tissues.

Roots are a good source to study the drought mechanism. To further confirm this, qRT-PCR showed the expression of TaRPK1 members in the leaves and roots of two-week-old seedlings with drought stress through PEG simulation. PEG-6000 treatment induced an upregulated expression in roots and leaf tissues in comparison to the susceptible genotype. Higher expression was observed in root seedlings in comparison to the leaf seedlings, except for TaRPK4 and TaRPK7, where higher expression was detected in the leaf tissues compared to the root tissues under drought stress (Figure 10). The TaRPK1 genes displayed higher expression in Pakistan 13 > Galaxy > Shafaq varieties. Furthermore, we also performed co-expression (Supplementary Figure S3) and interaction network (Supplementary Figure S4 and Supplementary Table S8) analyses and the results revealed that all the RPK1 members showed highly significant associations. These results indicate TaRPK1 gene involvement in drought stress regulation.