Full-length coding and protein sequences for the 16 CrRLK1L from rice were retrieved from the Rice Genome Annotation Project database (http://rice.plantbiology.msu.edu/). CrRLK1L homologs in wheat, hereafter denoted as Ta-CrRLK1L, were searched by using rice FL coding sequences in the NR database at NCBI (https://www.ncbi.nlm.nih.gov/) in Triticum aestivum. Sequences with high-scoring hits and maximum similarity with rice CrRLK1Ls were collected, and duplicates were eliminated. The NR sequences were used to search for respective A, B, and D homeologs in the Ensembl plants database (https://plants.ensembl.org/index.html), and genomic, CDS, and protein sequences for Ta-CrRLK1L with high-scoring hits and more than 95% identity were collected. Transcript variants that cover the full-length coding sequences were selected. Ta-CrRLK1L nucleotide coding sequences from Ensembl plants were used to search against the TSA and EST databases at NCBI for further sequence confirmation. The start codon position in the wheat genome for each Ta-CrRLK1L homeolog was determined by a blastn search of coding regions in the International Wheat Genome Sequencing Consortium (IWGSC) RefSeq v2.1 genome assembly at the WHEAT URGI database (https://wheat-urgi.versailles.inra.fr/). The physical properties, like the molecular weight of the protein and isoelectric point (pI), were calculated using the Expasy Compute pI/MW tool (https://web.expasy.org/compute_pi/). Subcellular localization for the proteins was predicted by WoLF PSORT (https://wolfpsort.hgc.jp/).

Protein sequences for genes encoding CrRLK1L for closely related species like Aegilops tauschii and Hordeum vulgare were retrieved from the NCBI database using balstp search with wheat protein sequences as queries. Protein sequences for Brachypodium distachyon were collected from Ensembl plants using T. aestivum nucleotide coding sequences as query, while sequences for Arabidopsis CrRLK1L gene families were retrieved from the TAIR database (https://www.arabidopsis.org/), and rice sequences were retrieved from the Rice Genome Annotation Project database (http://rice.plantbiology.msu.edu/). Sequences that share the highest sequence identity with wheat CrRLK1Ls were selected. E value cut-off of zero was used.

Conserved domains for the Ta-CrRLK1L proteins were confirmed by the Simple Modular Architecture Research Tool (SMART) (http://smart.embl-heidelberg.de/) and Batch Conserved Domain (CD)-Search tool (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) at NCBI, and active sites in the Ta-CrRLK1L domains were predicted by PROISTE (http://prosite.expasy.org/). Multiple sequence alignment for Ta-CrRLK1Ls was carried out by ClusatlW, and conserved domain features and active sites were visualized using the Jalview program (https://www.jalview.org/). For phylogenetic analysis, the multiple sequence alignment for protein sequences was carried out by ClustalW (https://www.ebi.ac.uk/Tools/msa/clustalo/), and the evolutionary relationships for CrRLK1L in T. aestivum and other five species, namely Ae. tauschii, A. thaliana, B. distachyon, H. vulgare, and O. sativa, were inferred by analyzing full-length amino acid sequences in MEGA11 (Tamura et al., 2021). Ta-CrRLK1L A homeolog copies were used in the analysis, except for Ta-CrRLK1L7, for which D copy was used. In total, eighty-eight CrRLK1L protein sequences were used in the analysis. Sequences were aligned by ClustalW with default parameters, and a phylogenetic tree was constructed using the Jukes-Cantor model by the Neighbour-joining method and 1000 bootstrap iterations. The phylogenetic tree was visualized by Tree of Life (iTOL) (https://itol.embl.de/login.cgi).

Gene duplication analysis for the Ta-CrRLK1L gene families was conducted using cut-off criteria that the sequence alignment had >80% coverage of that of the longer gene, and the aligned region had the identity of >80% at the nucleotide level among the gene pairs (Wang et al., 2016). The Ka/Ks (nonsynonymous substitution rate/synonymous substitution rate) ratio, which determines the selection pressure on the duplicated genes, was calculated for the paralogous gene pairs using TBtools (Chen et al., 2020). Gene duplications were illustrated using circos using shinyCircos-V2.0 (https://venyao.xyz/shinycircos/) (Wang et al., 2023). Chromosome length for each copy in wheat was collected from the Ensembl Plants database, and start and end positions for Ta-CrRLK1L genes were determined from the (IWGSC) RefSeq v2.1 genome assembly. The position of the Ta-CrRLK1L genes on the chromosomes was represented using MG2C v2.1 (http://mg2c.iask.in/mg2c_v2.1/) (Chao et al., 2021). Intron/exon junctions were determined by comparing Ta-CrRLK1L cDNA sequences from the NR database at NCBI with genomic sequences from Ensembl by using Splign (https://www.ncbi.nlm.nih.gov/sutils/splign/splign.cgi?textpage=online&level=form), and the intron/exon junctions were visualized by Gene Structure Display Server 2.0 (GSDS) (http://gsds.gao-lab.org/). Conserved motifs in CrRLK1L proteins were predicted using Multiple Em for Motif Elicitation (MEME) Suite 5.5.2 (http://meme-suite.org/tools/meme) with ten motifs as the number of motifs parameter and the GenomeNet Database (https://www.genome.jp/tools/motif/).

The nucleotide sequences from the upstream 2 Kb regions of the start codons of Ta-CrRLK1L genes were retrieved from the Ensembl Plant database and analyzed for cis-acting elements. Cis-acting elements were determined by using PlantCARE (https://bioinformatics.psb.ugent.be/webtools/plantcare/html/).

Tissue-specific gene expression and gene expression in response to stress conditions for Ta-CrRLK1L gene families were analyzed using T. aestivum transcriptome datasets available at SRA (Sequence Read Achieve) repositories in the NCBI database (https://www.ncbi.nlm.nih.gov/sra). To determine the tissue-specific relative gene expression levels, the transcriptome dataset of seventy-one tissues of wheat cultivar Azhurnaya, available at eFP Browser (http://bar.utoronto.ca/efp_wheat/cgi-bin/efpWeb.cgi), was used (Ramírez-González et al., 2018). Additionally, the gene expression levels in the different tissue types were analyzed from the tissue-specific expression data from BCS cv-1 Development in Wheatomics 1.0 database (http://wheatomics.sdau.edu.cn/expression/wheat.html).The relative gene expression levels in five tissue types, namely, inflorescence, leaf at the whole plant seed formation stage, root at the cotyledon emergence stage, seed at the fruit ripening stage, and stem, were analyzed from the transcriptome datasets available at SRA database at NCBI (Pingault et al., 2015). Tissue-specific expression in the stamen and pistil tissues was analyzed by using BioProject: PRJEB36244 datasets at SRA repositories (https://www.ncbi.nlm.nih.gov/sra). Differential gene expression in response to abiotic factors such as drought, heat, combined stress (Liu et al., 2015), cold (Li et al., 2015b), and salt stress induced by NaCl (Bio Project Accession: PRJNA632706) was carried out by using the transcriptome datasets available at SRA repositories at NCBI. The gene expression changes in response to infections by Claviceps purpurea, which causes ergot disease, and Fusarium graminiearum, responsible for head blight in wheat, were analyzed using transcriptome datasets obtained from SRA ( Supplementary Table S1 ).

The relative gene expression levels across five different tissue types and changes in gene expression in response to abiotic and biotic stresses were calculated as FPKM (Fragments Per Kilobase Per Million). The FPKM changes in gene expression were determined by search of Ta-CrRLK1L gene Ensembl Plant database identifiers against the T. aestivum RNA-seq Database (http://ipf.sustech.edu.cn/pub/wheatrna/) at the Plant Public RNA-seq Database (Zhang et al., 2020b), using restricted SRA accessions search for biological replicates for each dataset for five tissue types or abiotic or biotic stresses. This database is comprised of a pipeline that uses tools such as HISAT2 (Kim et al., 2015; Kim et al., 2019) to align raw reads to the reference genome and StringTie (Pertea et al., 2015) to assemble and quantify transcripts. The heatmap for FPKM values for the relative level of gene expression and the fold change in response to stress conditions were visualized by TBtools (Chen et al., 2020).

Gene co-expression network analysis for the moderately or highly expressed CrRLK1Ls in tissues was analyzed using tissue-specific RNA-Seq network in the Wheat co-expression network (WheatCENet) database using default parameters and positive co-expression. This database comprises of RNA-seq datasets for tissue specific gene expression in wheat (http://bioinformatics.cau.edu.cn/WheatCENet) (Li et al., 2023). The AgriGO v2.0 database (http://bioinformatics.cau.edu.cn/WheatCENet/GOanalysis.php) was utilized for the GO enrichment analysis of the co-expressed network nodes identified from the WheatCENet database using the Singular Enrichment Analysis (SEA) tool in the T. aestivum (Chinese spring species) with the default parameters and a False Discovery Rate (FDR) value of 0.05. The query consisted of the co-expressed network nodes Ensembl IDs of the genes.

In total, forty-three genes encoding CrRLK1Ls were identified in the hexaploid genome of Triticum aestivum, which consists of fifteen paralogous genes with three homeologous copies each, namely A, B, and D, from their progenitor species, except for Ta-CrRLK1L2 and Ta-CRLK1L7, which had missing -2-D and -7-A copies. However, the blastn search for Ta-CrRLK1L2 at IWGSC RefSeq v2.1 and Ensembl plants had a hit with 88% identity within the coding regions of the genes on Chr3A. To confirm the presence of a -3-D copy in the diploid progenitor, a blastn search in the Ae. tauschii genome in the NCBI database was done, which also did not show a high sequence similarity hit and suggested that -2-D may also be missing in Ae. tauschii. Ta-CrRLK1L2-A and 2-B had 789 aa and 814 aa length sequences in the Ensembl database, which were more like transcript variants, and these sequences were corrected using TSA and wheat WGA_v0.4 scaffolds at the IWGSC database. Among CrRLK1Ls, -7-B was the shortest paralog with a length of 630 aa and showed missing amino acids in the N-terminal region, which was also confirmed using the IWGSC database and whole genome shotgun (WGS) contigs at the NCBI database. Sequences for CrRLK1Ls were confirmed from at least two independent databases, including Ensembl Plants and the NR, TSA, and EST databases at NCBI. The transcript variants for Ta-CrRLK1L gene family members, -11-B (two variants), -4-B (two variants), and -14-D (4 variants), were identified from the Ensembl plant database, and the correct variant that codes for FL protein was selected for further analysis ( Table 1 , Supplementary Table S2A ).

The protein lengths and molecular weights for Ta-CrRLK1Ls were 801–911 aa and 88.7-99.2 KDa, respectively. Ta-CrRLK1L6-A and -13-B were translocated from 4A to 5A and 7B to 4A chromosomes, respectively. Most homoeologous copies of Ta-CrRLK1L showed more than 93.5% similarity in the protein-coding regions within their homeologs, except for -5-B and -5-D, which had 91.0% similarity, respectively. The sub-cellular localization predicted by using the WoLF PSORT tool showed that 36 Ta-CrRLK1Ls were localized to the plasma membrane. Surprisingly, -14-B, -15-A, -15-B, and -15-D were found localized to the vacuole, while -9-D, -12-A, and -12-D were localized to the chloroplast ( Table 1 ).

Blast searches using T. aestivum CrRLK1L protein-coding sequences in monocot species Ae. tauschii, B. distachyon, and H. vulgare identified 14, 11, and 15 CrRLK1L encoding genes, respectively. The CrRLK1L2-D copy was missing in Ae. tauschii ( Supplementary Table S2B ).

Multiple sequence alignment and conserved domain prediction using SMART and NCBI Batch-CD Search tool showed that T. aestivum CrRLK1Ls had highly conserved characteristic domains that included a malectin-like domain (PF12819) and a catalytic domain for serine/threonine protein kinases (SM000220) ( Figure 1A ). The lengths of the malectin-like domain in Ta-CrRLK1L varied between 270 and 378 aa, except for the -7-B gene copy, which had a 198 aa malectin-like domain. This suggests that the point mutation caused the truncated malectin-like domain in the -7B copy of Ta-CrRLK1L. The serine/threonine protein kinase domain lengths were ranged from 265 to 281 aa. Sequence analysis using PROSITE showed that Ta-CrRLK1Ls were comprised of featured domains that include an ATP-binding region signature (PS00107) and an active site signature for serine/threonine protein kinases (PS00108). The serine/threonine protein kinases had the active site residue D, which is directly involved in the catalytic functions predicted by PROSITE ( Figure 1B ).

Phylogenetic analysis of CrRLK1Ls using full-length protein sequences showed that more closely related species like T. aestivum, Ae. tauschii, and H. vulgare were clustered together. The paralogs of Ta-CrRLK1Ls were clustered with Arabidopsis CrRLK1Ls, which have been known to be involved in plant reproduction and pollen-pistil interaction in Arabidopsis. At-FER was clustered with Ta-CrRLK1L5, -9, At-ANX1 and At-ANX2 were clustered with -12, At-HERK2 was clustered with -10, At-BUPS1 and At-BUPS2 were clustered with -13, At-CAP was clustered with -1, and At-HERK1 and At-ANJEA were clustered with -4, -8, -11 and -14 ( Figure 2 , Supplementary Figure S1 ), which suggests that these genes may also have a similar role in plant reproduction.

Phylogenetic analysis in T. aestivum revealed clustering of closely related CrRLK1Ls such as -1 and -13, -5 and -9, -6 and -7, -10 and -15, as well as -11 and -14 ( Figure 3A ). The intron/exon junctions determined by Splign showed that most Ta-CrRLK1Ls had no introns and were comprised of single exons, except for the paralog -11-D and homeologs of -14, which have a single intron and two exons, suggesting the structural diversity within gene family ( Figure 3 , Supplementary Table S3A ). Ta-CrRLK1-11-D copy had the smallest and largest introns of 124 bp and 1149 bp, respectively. The homeolog of -14 had shorter first exons of 36 bp, 36 bp, and 111 bp for -14-A, -14-B, and -14-D, respectively. The MEME suite predominantly predicted motifs associated with serine/threonine protein kinase domains ( Figures 3C, D , Supplementary Table S3B ).

Chromosome 4A, 4B, and 4D had the highest distribution of CrRLK1L genes and showed the distribution of four genes, including homeologs for Ta-CrRLK1L5, -8 and -12, and paralogs -6B and -6D. However, 4A showed the translocated -13-B paralog from 7B to 4A chromosome ( Figure 4A , Table 1 ). Chromosomes 3D and 7B had a distribution of single CrRLK1L genes that included -4-D and -14-B, respectively. Chromosome 6 did not have any distribution of CrRLK1L genes and, hence, was eliminated ( Figure 4A ).

Gene duplication analysis using the criteria of an 80% cut-off for alignment coverage and sequence identity at nucleotide level among Ta-CrRLK1L gene families showed that only the gene pair CrRLK1L6 and -7, and -5 and -9 paralogs had the gene duplication events ( Figure 4B ). Other CrRLK1L paralogs did not meet these criteria. The Ka/Ks ratio analysis suggested that these duplicated gene pairs had a value of less than 0.45, indicating that the evolution of these CrRLK1L genes is under purifying selection ( Supplementary Table S4 ).

Gene expression analysis for CrRLK1Ls in T. aestivum across seventy-one tissues in Azhurnaya spring wheat showed that CrRLK1L-12 and -13 had higher expression levels in anther tissues. The homeolog -13-D exhibited the highest expression level at 158.4 TPM, making it the most highly expressed paralog in the CrRLK1L gene family. Similarly, paralog -12-B also showed a higher expression level of 61 TPM in anther tissues. ( Figure 5A , Supplementary Table S5A ). Gene expression for different tissue types analyzed from the Wheatomics database (http://wheatomics.sdau.edu.cn/) were also in agreement with the gene expression levels of -12 and -13 in anther as well as other tissues ( Supplementary Table 5B ). In the phylogenetic analysis, Ta-CrRLK1L-12 grouped closely with At-ANX1 and At-ANX2, whereas Ta-CrRLK1L-13 clustered with At-BUPS1 and At-BUPS2 ( Figure 2 , Supplementary Figure S1 ). These genes have a role in maintaining pollen tube integrity and preventing premature rupture before reaching synergid cells. This suggests that Ta-CrRLK1L-12 and -13 may have similar functions in wheat. Ta-CrRLK1L5, -9, and -14 exhibited moderate expression levels in the stigma and ovary tissues and were clustered with At-FER, whereas -11 and -14 were clustered with At-HERK1 and At-ANJ, known for their roles in pollen tube reception at synergid cells during reproduction, suggesting their involvement in pollen tube reception during fertilization.

The relative level of gene expression analyzed across five tissue types was in agreement with the analyzed expression for seventy-one tissue types for CrRLK1L3 and showed that it was highly expressed in leaf at whole plant fruit formation (30–50% moisture stage), with expression values of 29.2 to 36.3 FPKM ( Figure 5B , Supplementary Table S6 ). Ta-CrRLK1L5 and -14 were expressed across five tissue types with different expression levels. Paralog -5-B was least expressed in seed with 8.8 FPKM, while -5-D was highly expressed in stem with 34.5 FPKM. Homeologs of -6 were expressed in stem, leaf, and root tissues from 3.43 to 13.89 FPKM.

To explore Ta-CrRLK1L’s role in plant reproduction, gene expression was analyzed in immature stamen and pistil tissues of a cultivar Fielder. Ta-CrRLK1L5, -9, and -14 exhibited higher expression than other family members in both stamen and pistil tissues; however, the gene expression levels in pistil tissues were higher than in stamen tissues. -5-A and -14-D had a higher level of expression than their homeologs, with 38.3 and 29 FPKM in pistil tissues, respectively, which was in agreement with the seventy-one tissues transcriptome data ( Figure 5C , Supplementary Table S7 ). Similar paralogs were also expressed in stamen tissues at 27 and 31 FPKM values, suggesting dual roles of these genes in male as well as female side in pollen-pistil interaction.

Gene co-expression network analysis of seven moderately or highly expressed paralogs, namely -3A, -5A, -9A, -11A, -12A, -13A, and -14A using the WheatNet database revealed the functional and regulatory roles in morphogenesis and developmental processes ( Supplementary Tables S8A–C ). GO enrichment analysis for -3A co-expression network nodes showed the enrichment in enzyme linked receptor protein signaling pathway (GO:0007167), plant epidermal cell differentiation (GO:0090627), root epidermal cell differentiation (GO:0010053) and cell surface receptor signaling pathway (GO:0007166). Co-expression nodes for -5A showed the enrichment mostly in the signaling pathways such as cellular response to organic cyclic compound (GO:0071407), transmembrane receptor protein tyrosine kinase signaling pathway (GO:0007169), and steroid hormone mediated signaling pathway (GO:0043401). Co-expression network of -9A showed enrichment in cellular component organization (GO:0016043), embryo sac development (GO:0009553), ubiquitin protein ligase binding (GO:0031667), symplast (GO:0055044) and External encapsulating structure (GO:0030312). Similarly, for -12A and -13A co-expression network enrichment in the developmental growth (GO:0048589), unidimensional cell growth (GO:0009826), cell growth (GO:0016049) and developmental growth involved in morphogenesis (GO:0060560O) was observed, suggesting the redundancy in the functioning of these CrRLK1Ls. These results suggested that these CrRLK1Ls control the cell morphogenesis and developmental processes through different pathways while others act redundantly though similar pathways. The details of the positive interactions of the CrRLK1Ls with the co-expressed genes, and co-expression nodes, edges, and GO analysis is given ( Supplementary Figure 2 and Supplementary Tables S8A-C ).

Most of the CrRLK1L paralogs were upregulated or downregulated in response to one or more abiotic factors. In response to the cold treatments, the paralogs Ta-CrRLK1L4-B and -15-A were up-regulated by 3.0 and 3.3-fold, respectively. The paralogs that were downregulated to 15%, 18%, and 21% compared to control conditions included -3-A, -7-B, and -11-A. CrRLK1L-3-A had the highest 18% decrease in gene expression, from 24.69 FPKM to 4.42 FPKM. Other paralogs had a 30–47% decrease in the gene expression level in response to cold stress, which included one of the homeologs of Ta-CrRLK1L3, -4, -6, -7, and -11 ( Figure 6A , Supplementary Table S9 ).

The Ta-CrRLK1L gene family members showed changes in gene expression levels in response to drought, heat, and combined stress. The homeologs of -9 and paralogs -10-A, and -14-B were induced by more than 2.3-folds in response to drought stress of 1 hr and 6 hr, while homeologs of -2, -4, -7, -9, -10, and -15 were downregulated in the same, with 4-B exhibiting highest downregulation of 25% in the 6 hr drought stress treatment. The paralogs upregulated in response to 1 hr or 6 hr of heat stress included -4-B, -9-A, -10-D, -11, -15-B and -15-D. In contrast, -2, -3, -4, 6, -7, -9, -10, and -11 were downregulated, with the highest downregulation for 10-B from 1.04 to 0.04 FPKM. In combined stress, -15-B, was induced 13.9-fold, while 10-B was downregulated from 1.04 FPKM to 0.02 FPKM. The homeologs -7-B and -7-D were downregulated to 12% of their control levels, with FPKM values of 7.28 and 3.94 ( Figure 6B , Supplementary Table S10 ).

Changes in gene expression analyzed in response to 100 mM NaCl treatments Ta-CrRLK1L-1-A had the highest induction of 5.77-fold compared to -1-B and -1-D, which were induced by 3.35 and 2.93-fold, respectively. However, the lower expression levels of 0.12, 0.08, and 0.23 FPKM for these homeologs were detected in control treatments. These results suggest that Ta-CrRLK1L1 may have a role in salt stress responses ( Figure 6C , Supplementary Table S11 ).

In response to Fusarium graminearum (Fg) inoculation, Ta-CrRLK1L-7-D was highly induced from 0.08 FPKM in control to 0.81, 1.85, and 1.75 FPKM in 24 hr, 48 hr, and 72 hr post-inoculation treatments. -6-D was induced by a 6.93-fold increase in 48 hr post-Fg inoculation. However, in response to 72 hr post-inoculation, -15 homeologs were downregulated to 10-12% of their control values ( Figure 6D , Supplementary Table S12 ).

Claviceps purpurea (Cp) infection in the stigma induced -6-D and -15-D paralogs from 0.06 FPKM in mock to 0.50 and 0.57 FPKM in 1 hr post-inoculation treatment, whereas the homeologs -12-A and -12-D were downregulated to 14% and 11% in the treatment compared to their control ( Figure 6E , Supplementary Table S13 ). The paralogs -3-A and -6-D were induced by 10.26 and 10.92-fold in response to 48 hr post Cp inoculation in the ovule. However, in the control conditions, -3-A and -6-D have lower expression of 0.14 and 0.06 FPKM, respectively. -6-D and -3-A were also upregulated 48 hr post-inoculation in the ovary and transmitting tissue to a fold change of 7.04 and 4.26, respectively.

The prediction of cis-acting elements in the putative promoter regions of CrRLK1L genes revealed that 54% consisted of core promoter and enhancer region elements, including TATA and CAAT boxes. The remaining 46% was categorized into three groups, namely, cis-acting elements associated with hormonal responses, plant developmental processes, and abiotic and biotic stress conditions ( Figure 7A ). These categories were further subcategorized. For instance, hormonal responses were divided into abscisic acid, auxin, gibberellin, and salicylic acid ( Figure 7B ), while other elements were grouped based on their functions ( Figures 7C, D ). All Ta-CrRLK1L genes are comprised of cis-acting elements for light and abscisic acid. Ta-CrRLK1L1 had cis-acting elements for palisade mesophyll cell differentiation, and -14 had meristem-specific activation elements, which fell under the plant developmental processes category. The paralogs -9-B, -12-A, -12-D, -13-A, and -13-B possessed elements for seed-specific regulation. Cis-acting elements associated with defense and stress were predicted in all CrRLK1L genes except -10, -11, and -14. Additionally, at least one of the homeolog CrRLK1L genes showed cis-acting elements for methyl jasmonate ( Supplementary Table S14 ).