To identify LRR-RLK genes in T. aestivum, the amino acid sequences of Arabidopsis (AtLRR-RLKs) and rice LRR-RLKs were used as queries to perform BLAST searches in Ensembl Plants database (http://plants.ensembl.org/) (Bolser et al., 2015). Additionally, we downloaded the corresponding hidden Markov models (HMM) of LRR and PK domains from the Pfam database version 32.0 (http://pfam.xfam.org/) (Finn et al., 2014) and used them as references to search for LRR-RLKs in the T. aestivum protein sequences from the International Wheat Genome Sequencing Consortium (release-55) (http://plants.ensembl.org/index.html). The results of these two approaches were combined to identify the TaLRR-RLKs wheat genome with an E-value < 10^-5. All candidates were further screened via eliminating the redundant sequences and functional annotations and then analysis with the NCBI Batch Web CD-Search Tool (https://www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi). The presence of signal peptides (SP) was calculated with SignalP (http://www.cbs.dtu.dk/services/SignalP/) (Petersen et al., 2011), the prediction of transmembrane (TM) domains was performed using the TMHMM server (http://www.cbs.dtu.dk/services/TMHMM/) (Krogh et al., 2001) and the subcellular localization of TaLRR-RLKs were predicted using ProtComp software on the SoftBerry online database (http://linux1.softberry.com/) (Long et al., 2014), respectively.

Multiple sequence alignment of AtLRR-RLKs and TaLRR-RLKs was carried out using MEGA6 software (Tamura et al., 2013). Furthermore, an unrooted phylogenetic tree was generated with the maximum likelihood (ML) method by 1 000 bootstrap tests (Price et al., 2010). The TaLRR-RLKs were classified according to the topology of the phylogenetic trees and the established classification rely on Arabidopsis homologs (Shiu and Bleecker, 2001b). Pairwise deletion was used for data processing, while the Poisson distribution was adopted for tree-building model. The online tool iTol (http://itol.embl.de) (Letunic and Bork, 2011) was used to create improved graphical representations of the phylogenetic tree.

To investigate the distribution of TaLRR-RLKs in wheat chromosomes and gene duplication events, we downloaded wheat genomic sequences and genome annotation files from the Ensembl Plants database for further analysis. Subsequently, the chromosomal locations and synteny relationships of TaLRR-RLKs were illustrated using TBtools software using the Graphics function (Chen et al., 2020).

To better understand the protein features of each subfamily, all the identified TaLRR-RLKs were searched in Pfam database (Finn et al., 2014), and verified using CDD analysis (http://www.ncbi.nlm.nih.gov). To characterize the structural divergence of TaLRR-RLKs in each group, the conserved motifs in their encoded proteins were identified using TBtools software and the Simple Multiple Em for Motif Elicitation (MEME) Wrapper function (Chen et al., 2020). The maximum number of motifs was set to 10 and other optional parameters to default. The results of conserved domains and motif analysis of several representative proteins from each subfamily were displayed using TBtools (Chen et al., 2020), and then mapped in the clustered phylogenetic tree using Photoshop software. We used the automated SWISS-MODEL homology modeling server (https://swissmodel.expasy.org/) to generate three-dimensional protein models of the TaLRR-RLK proteins (Waterhouse et al., 2018). We randomly picked a gene from each group to conduct homology modeling based on QMEAN (-4 to 0, where the closer the value to 0, the better the match between the examined protein and the template) (Benkert et al., 2011) and GMQE (0 to 1, where the larger the value, the better the quality) (Waterhouse et al., 2018).

The 2000 bp upstream of all TaLRR-RLKs coding sequences were identified from the Ensembl Plants database and then used predicted the conserved cis-acting elements in the putative promoter regions by the PlantCARE software (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) (Lescot et al., 2002). All identified cis-acting elements were clustered and arranged according to their different functions and visualized with the TBtools software (Chen et al., 2020).

T. aestivum. cv. Yangmai 158 was grown in a glasshouse at 15 ± 2°C with a photoperiod of 16 h light, 8 h dark and 70% relative humidity (Yang et al., 2020). To analyze the expression profile of TaLRR-RLKs in wheat, we extracted total RNA from seven tissues (from bottom to top: root, stem, first leaf, second leaf, third leaf, fourth leaf, and fifth leaf) of five-leaf stage wheat seedlings for quantitative real-time PCR (qRT-PCR). To analyze the expression profile of TaLRR-RLKs under hormone treatments, wheat seedlings with two-leaf-stage were sprayed with 100 μmol L^–1 abscisic acid (ABA), 100 μmol L^–1 methyl jasmonate (MeJA), or distilled water (as control) (Yu et al., 2019). Three biological replicates of samples were collected at four different time after hormone treatment (2, 4, 6 and 12 h) for total RNA extraction. For investigate the expression level of TaLRR-RLKs in wheat leaves, the full-length cDNA clones of CWMV were used for CWMV inoculation according to a previous study (Yang et al., 2016). At 7 days post inoculation (dpi), the systemic leaves were sampled for RNA extraction and qRT-PCR assay.

The VIGS tool function on the online website (https://solgenomics.net/)was used to obtain the best VIGS fragment sequence (300 bp) of TaLRR-RLKs. The specific fragments were amplified from wheat cDNA using RT-PCR ( Supplementary Table S1 ), digested with PacI and NotI restriction enzymes, and inserted into the BSMVγ gene (BSMV : TaLRR-RLK). Plasmids BSMV α and γ were linearized with MluI restriction enzyme, BSMV β with SpeI, and BSMV γ-gene with BssHII. Then, the linearized plasmid DNAs of BSMV RNAα, β, and γ were prepared for transcription in vitro using the T7 in vitro transcription kit (Ambion, Austin, TX; Promega, Shenzhen, China). Next, the transcription products were mixed at a 1:1:1 molar ratio, followed by addition of excess FES buffer (0.06 M potassium phosphate, 0.1 M glycine, 1% bentonite, 1% sodium pyrophosphate decahydrate, 1% celite, pH 8.5) for friction inoculation into second leaves of two-leaf-stage wheat plants, as previously described (Yang et al., 2020). Wheat seedlings inoculated with FES buffer were used as negative controls. Finally, the inoculated wheat was placed in a growth environment at 28 °C and 70% relative humidity for 7-10 days to observe the symptoms. Similarly, we mixed the transcripts of linearized plasmid DNAs of CWMV RNA R1 and R2 at a 1:1 molar ratio and inoculated third leaves of the BSMV-infected wheat plants. Subsequently, the virus inoculated wheat plants were grow for 21 days at 15 ± 2 °C and 70% relative humidity. At 10 dpi, viral infection leaves were harvested for RNA extraction. Photographs of viral infection leaves were taken at 40 dpi.

Total RNA was extracted from each sample using a HiPure Plant RNA Mini Kit (Magen, Guangzhou, China). First Strand cDNA Synthesis Kit (Toyobo, Kita-ku, Osaka, Japan) with random primers was used to synthesize first-strand cDNA, and 1 µg total RNA was added per 20 µL reaction volume, as previously described (He et al., 2020). Hieff qPCR SYBR Green Master Mix (Yeasen, Shanghai, China) and ABI 7900HT sequence detection system (Applied Biosystems QuantStudio 5, Foster City, CA, USA) were used to perform quantitative PCR. The 2^-ΔΔCt method (Livak and Schmittgen, 2001) was used to quantify the relative gene expression. In addition, the expression of T. aestivum cell division cycle (TaCDC) gene (Accession Number: XM_020313450) was analyzed and used as an internal control (Zhang et al., 2019). Each treatment in this study had at least three biological and three technical replicates. All the gene-specific primers used in this study were listed in Supplementary Table S1 .

Microsoft Excel was used to determine the mean values and standard errors of the treatments. t-test or the Tukey’s test was performed using the SPSS 16.0 software (SPSS, Inc., Chicago, IL) to determine the significance of differences. The significant difference with unequal variance between two treatments was decided using the probability (P) value <0.05.

For purpose of identify the LRR-RLK family members in wheat, we used LRR-RLK sequences from the A. thaliana and rice genomes as a query to identify LRR-RLK sequences in the wheat genome database using the BLASTP and HMM approach. Total 6967 identified proteins were then submitted to the Pfam and SMART databases for analysis of domain structures. 929 candidates with at least two LRR domains and a kinase domain in the wheat genome were identified. The detailed information, including the accession number of each protein, were listed in Supplementary Table S2 . To study the evolutionary relationships of TaLRR-RLKs, all 929 candidates and all AtLRR-RLK proteins was used to construct an ML tree ( Figure 1 ). The results of the phylogenetic analysis illustrated that all TaLRR-RLKs were clustered in 14 major groups according to the nomenclature of the Arabidopsis homologues within the same group (I-XI, XII-a, XII-b, XIII) ( Figure 1 ; Supplementary Table S2 ). Of these, subfamilies XII-a, III, XI, and XII-b were the top four largest subfamilies, with 305, 109, 96, and 93 members, respectively. The subfamilies IV and IX were the two smallest, with 14 and 7 members, respectively. Similar to phylogenetic tree in Arabidopsis, the members in LRR-XII fell into two different subfamilies (LRR-XIIa and LRR-XIIb). Although a few members belong to subfamily XII among AtLRR-RLKs, total 398 TaLRR-RLKs were identified as orthologs of these AtLRR-RLKs ( Figure 1 ).

Based on the initial position of each gene on the wheat chromosomes, 929 TaLRR-RLK genes were randomly distributed among all 21 chromosomes (chromosomes 1A to 7D, Chr1A to Chr7D), with Chr2D having the most (65) and Chr4B having the least (30) genes ( Figure 2 ; Supplementary Table S3 ). There were 16 TaLRR-RLKs identified on the unknown wheat chromosome (ChrUn). Combined with wheat gene density, the distribution of TaLRR-RLKs was not significantly correlated with the distribution of wheat genes. To explore duplication relationship of TaLRR-RLKs during evolution, we performed synteny analysis according to the position of each gene in wheat chromosomes. In our analysis, 123 genes were tandemly duplicated ( Supplementary Table S4 ), and a total of 921 collinear events were identified and distinguished in different colors ( Figure 2 ; Supplementary Table S5 ).

To gain insight into the spatial conformations of TaLRR-RLKs, we used the SWISS-MODEL template library to conduct homology modeling and generate three-dimensional protein models. We randomly selected a TaLRR-RLK from each subfamily as a representative to analyze its predicted structure and illustrate the results in Figure 3 . All 14 TaLRR-RLKs could be forecast as models, revealing that they maintained structural integrity during evolution, which is important for their function. Meanwhile, we found that the spatial structures of proteins from the same clade are highly similar.

The deduced structures of representative proteins of each subfamily were mapped in the phylogenetic tree ( Figure 4A ), and the exact distributions of the LRR domains and PK domains identified for each TaLRR-RLK were listed in Supplementary Tables S6 , S7 , respectively. Our results showed that all members of each LRR-RLK subfamily contain multiple types of LRR domains. The frequency of LRR repeats was similar in each LRR-RLK subfamily and PK domains did not vary significantly across different subfamilies ( Figure 4B , Supplementary Tables S6 , S7 ). The most common types of LRR domain were LRR_8 (leucine-rich repeat 8), LRR_1 (leucine-rich repeat 1), and LRRNT_2 (leucine-rich repeat N-terminal 2), with different distribution patterns in each subfamily. For instance, subfamily IX preferentially contained LRRNT_2 but not LRR_1 and two LRR_8 domains were identified in this subfamily. Subfamily I preferentially contained LRR_8 and LRR_4 (leucine-rich repeat 4) but not LRRNT_2 or LRR_6 (leucine-rich repeat 6). Additionally, four subfamilies, including IV, VI, VIII, and IX did not contain LRR_4 or LRR_6. To further elucidate the structural characteristics and potential functions of LRR-RLKs in wheat, we set the threshold to 10 motifs in the MEME online tool to predict the putative motifs of these proteins. The distribution of these putative motifs in each representative TaLRR-RLK was shown in Figure 4C , and the detailed distribution was listed in Supplementary Table S8 . Among them, motifs M6, M5, M1, M3, M8, and M2 were located in the PK domain. Motifs M1, M2, M3, M4, M5, and M6 were present in all members of each subfamily. The M4 motif was frequently detected in all TaLRR-RLK subfamilies, particularly in subfamilies XI, XII-a, and XII-b. The M9 motif had a dramatic presence in subfamilies I, II, IV, and VII. M1 appeared most frequently in subfamily XII-b. Meanwhile, we submitted protein sequences of each TaLRR-RLK subfamily to the SignalP website to search for possible signal peptides and found that 706 of the 929 TaLRR-RLKs contained SP, whereas 39% of the subfamily XII members did not contain signal peptides ( Supplementary Table S9 ). Additionally, 865 genes contained at least one TM domain based on TMHMM prediction. Subcellular localization analysis showed that 900 TaLRR-RLKs were localized to the plasma membrane ( Supplementary Table S9 ).

To analysis of promoter elements of TaLRR-RLKs, we used PlantCARE to detected cis-acting regulatory elements of the promoter region (2000 bp sequence upstream of the translation start site) of the 929 LRR-RLKs. The results showed that a total of 80 putative cis-acting elements were prevalent in TaLRR-RLKs. The details of all genes was shown in Table 1 and Supplementary Table S10 . Briefly, the cis-acting elements predicted in the promoter regions were divided into seven categories ( Table 1 ). The most abundant elements were phytohormone regulation-related cis-acting elements that respond to auxin, ABA, MeJA, gibberellins (GA) and salicylic acid (SA). The MeJA response element mainly included the CGTCA- and TGACG-motifs, and the ABA responsive element consisted of ABRE. Moreover, the TGA-element, TGA-box, AuXRR-core and AuxRE were involved in auxin responsiveness and the GARE-motif, P-box, and TATC-box were considered gibberellin-responsive elements involved in the regulation of gibberellin responsiveness. The second type of widely distributed regulatory elements were light responsive elements, including Box 4, G-Box, GATA-motif, TCT-motif, TCCC-motif, GT1-motif, and Sp1, of which G-Box appeared to be the most widespread. Moreover, we also identified cis-acting regulatory elements related to environmental stimuli, including abiotic and biotic stresses, of which LTR and TC-rich repeats were the most abundant ( Table 1 ; Supplementary Table S10 ).

Based on the results of cis-acting regulatory element analysis, 12 types of hormone-responsive elements were identified ( Table 1 ; Supplementary Table S10 ). Among them, the hormone-responsive elements related to MeJA and ABA were widely distributed in the promotor region of TaLRR-RLKs ( Supplementary Table S10 ). To better understand the underlying functions of TaLRR-RLK response to MeJA and ABA, one TaLRR-RLK from each clade was randomly selected to analyze its respective expression level after hormone treatment using qRT-PCR. As shown in Figure 5 , TraesCS5D02G161800(I), TraesCS1A02G240300(III), TraesCS4B02G055600(IV), TraesCS2D02G034200(X), and TraesCS7D02G365800(XI) were insensitive to MeJA and ABA. By contrast, the expression level of TraesCS5B02G365500(V) and TraesCS4D02G172000(VII) was significantly induced by MeJA, while that of TraesCS1A02G391700(XII-a) and TraesCS4B02G211800(XIII) was reduced. The expression levels of TraesCS3B02G459200(VIII) were significantly up-regulated at 2 hours post treatment (hpt) and 6 hpt, but down-regulated at 4 dpt by MeJA. Besides, the expression levels of TraesCS2D02G321400(II) and TraesCS4D02G050200(IX) were up-regulated after treatment with ABA, while TraesCS5A02G428900(VI) and TraesCS1D02G002900(XII-b) were down-regulated. The qRT-PCR values used for the heat map are listed in Supplementary Table S11 . These results imply that TaLRR-RLK genes of T. aestivum participate in a variety of hormone responses ( Figure 5 ).

The plants were divided into seven different tissues: root, stem, first leaf, second leaf, third leaf, fourth leaf, and fifth leaf (bottom-top) to comprehensively explore the expression patterns of the TaLRR-RLKs. As shown in Figure 6 , taking the fifth leaf as a control, all 14 random selected genes from each subfamily were differentially expressed in each tissue and most tended to be down-regulated in the roots. The expression level of TraesCS5B02G365500(V) was down-regulated in each tissue. In contrast, the expression level of TraesCS4D02G050200(IX) was high in all tissues except for the root. In addition to up-regulation in the stem, TraesCS4B02G211800(XIII) showed moderate expression levels in other tissues. The expression levels of TraesCS1A02G240300(III) showed to be significantly up-regulated in the root and stem. TraesCS1D02G002900(XII-b) was observed to be up-regulated in the first, second, and third leaves, whereas its expression level in the fourth leaf was down-regulated in the root and stem ( Figure 6 ). The qRT-PCR values used for the heat map are listed in Supplementary Table S11 .

To investigate the function of TaLRR-RLKs in wheat resistance to CWMV infection, we analyzed the gene expression profile of TaLRR-RLKs under CWMV infection. The results indicated that four of selected 14 TaLRR-RLKs were differentially induced by CWMV infection, including TraesCS5D02G161800(I), TraesCS4D02G172000(VII), TraesCS7D02G365800(XI) and TraesCS4B02G211800(XIII). Besides, the accumulation level of TraesCS1D02G002900(XII-b) was significantly decreased under CWMV infection ( Figure 7 ). Based on the homologous analysis of these five genes, we renamed TraesCS4D02G172000(VII) as TaIRK, TraesCS7D02G365800(XI) as TaRLK5, TraesCS4B02G211800(XIII) as TaFEI1, TraesCS5D02G161800(I) as TaRLK-I.1 and TraesCS1D02G002900(XII-b) as TaRLK-XIIb.1, respectively. To further investigate the biological function of TaRLKs in CWMV infection, a BSMV-mediated gene silencing assay (BSMV-VIGS) was used to silence these five TaLRR-RLKs, respectively. We inoculated two-leaf-stage wheat seedlings with either BSMV:00+CWMV or BSMV : TaLRR-RLKs+CWMV. RT-PCR analysis verified the successful infection with BSMV and CWMV in all co-inoculated plants at 7 dpi (days post inoculation) ( Supplementary Figure S1 ). After 40 dpi, the mosaic symptoms were observed in leaves of all virus infected wheat plants, and the mosaic symptoms that appeared on BSMV : TaFEI1+CWMV and TaRLK-I.1+CWMV co-inoculated plants were significantly stronger than those on BSMV:00+CWMV co-inoculated plants. BSMV : TaIRK+CWMV co-inoculated plants exhibited much less mosaic symptoms and BSMV : TaRLK5+CWMV and BSMV : TaRLK-XIIb.1+CWMV co-inoculated plants showed similar mosaic symptoms to BSMV:00+CWMV co-inoculated plants ( Figure 8A ). Moreover, qRT-PCR analysis with specific primers confirmed that the five TaLRR-RLKs were successfully silenced in all BSMV : TaLRR-RLKs+CWMV co-inoculated plants ( Figure 8B ). Furthermore, the accumulation of CWMV CP was significantly up-regulated in TaFEI1-silencing plants, but down-regulated in TaIRK-silencing plants ( Figure 8C ). These results indicate that TaLRR-RLKs play an important role in wheat resistance to CWMV infection. Because CWMV RNA accumulation was significantly induced in TaFEI1-silencing plants compared to that in other selected LRR-RLKs silencing plants, indicating that TaFEI1 is more important for wheat resistance against to CWMV infection ( Figures 8A , C ). TaFEI1 was selected for further study.

Considering that the expression levels of TaFEI1 were significantly down-regulated after MeJA treatment, but not after ABA treatment ( Figure 5 ; Supplementary Table S11 ), we hypothesized that the function of TaFEI1 in CWMV infection were related to MeJA pathways. Thus, we calculated the accumulation of MeJA-response gene TaHOX24 and CWMV RNA in TaFEI1-silencing plants after MeJA hormones treatment using CWMV CP specific primers. The results of qRT-PCR showed that TaHOX24 was significantly induced in BSMV: 00 inoculated wheat leaves, but reduced in TaFEI1-silencing plants. However, MeJA treatment could not alter the relative expression level of TaHOX24 in TaFEI1-silencing plants ( Figure 8D ). TaFEI1-silencing plants was then inoculated with CWMV together with MeJA treatment. Compared to wheat BSMV inoculated plants without MeJA treatment, CWMV RNA was significantly reduced by MeJA in BSMV: 00-inoculated plants, but not changed under MeJA treatment in BSMV: TaFEI1-inoculated plants ( Figure 8E ). At 40 dpi, BSMV: 00+CWMV co-inoculated wheat leaves under MeJA exhibited much less symptoms of mosaic disease compared to that without MeJA. However, MeJA treatment was not altered the symptoms of mosaic disease TaFEI1-silencing plants ( Figure 8F ). These results suggested that TaFEI1 may be involved in wheat resistance to CWMV infection dependent on MeJA signal pathway.