Genome-wide association studies of Verticillium wilt resistance in a core collection of G. hirsutum

To identify loci responsible for VW resistance in upland cotton, we collected the VW resistance phenotypes from a core collection of G. hirsutum germplasm resources in the disease nursery at two different environments Anyang and Shihezi. Both DP and DI were used to evaluate the VW resistance of each accession in the panel. The terms AYDP and SHZDP represent the DP in Anyang and Shihezi, respectively. Similarly, AYDI and SHZDI represent DI in Anyang and Shihezi. The Pearson correlation showed strong correlation (r values ranging from 0.59 to 0.88) among the four traits (Additional file 1: Fig. S1c).

We next performed GWAS on AYDP, AYDI, SHZDP, and SHZDI using EMMAX [27] to identify the loci underlying each trait. The latest high-quality G. hirsutum TM-1 genome was used as the reference genome (CRI-TM-1, version1.0) [28]. First, all the traits’ genome-wide significance thresholds were set to a uniform threshold (5.03E − 6, P = 1/n, where n is the effective number of independent SNPs). As a result, 4 loci (AYDP1, AYDP2, AYDP3, and AYDP4) from A10, A13, D03, and D07, respectively, were identified for AYDP. For SHZDP, 2 loci (SHZDP1 and SHZDP2) from A08 and A10, respectively, were identified. For AYDI and SHZDI, 3 and 5 loci were detected respectively (Fig. 1a, Additional file 1: Fig. S2b and S3a, and Additional file 2: Table S1). Remarkably, AYDP1, AYDP2, and AYDP3 overlapped with previously reported QTL regions associated with VW or Fusarium wilt resistance (Fig. 1a) [15, 29]. Additionally, SHZDI1 was co-located with SHZDP2 and AYDP1, SHZDI3 was co-located with AYDP4, and AYDI3 was co-located with AYDP2, which provides important loci for targeted-genotyping-based selection.

Fig. 1

Genome-wide association using SNPs and introgression analysis on SHZDI1/SHZDP2/AYDP1. a Manhattan plot for GWAS on the Verticillium wilt for disease percentage (DP) and disease index (DI) in samples from Shihezi and Anyang. High-confidence QTLs and previously reported QTLs are marked in black and red words, respectively. b Fst analysis of A10 between the tolerant group and susceptible group based on lead SNP of SHZDI1. c SNP schematic around SHZDIP1 showing 2 haplotypes (Hap. 1 and Hap. 2) in G. arboreum and 5 haplotypes (Hap. A to Hap. E) in G. hirsutum. Number of materials is marked on the left. Hap.B of G. hirsutum carrying the introgressed G. arboreum fragments of Hap.1. Blue box indicates an introgressed region. d DI of haplotypes in greenhouse-grown G. arboreum and G. hirsutum. Significant difference was measured between different groups (Student’s t test). G. hirsutum carrying the introgressed fragments showed significantly enhanced VW resistance compared to materials without these fragments. G. arboreum variants carrying either Hap. 1 or Hap. 2 were more tolerant than G. hirsutum. e SNP-based phylogenetic tree showing that Hap.B lines of G. hirsutum harbor the introgression from A10 of G. arboreum

SHZDI1 originated from an exotic introgressed fragment of G. arboreum

Based on the allelic genotypes of the lead SNP of SHZDI1/SHZDP2/AYDP1 (Additional file 1: Fig. S1d), we divided the GWAS panel into 2 groups: tolerant and susceptible. We found that the Fst value between the two groups significantly increased on A10 ranging from 112.54 to 114.07 Mb (Fst peak = 0.982) (Fig. 1b). This indicates that the tolerant and susceptible groups demonstrated strong differentiation in this region. Since exotic introgression can result in strong differentiation, we speculated whether this observation is caused by introgression or not.

The triple hybrid (G. arboreum × G. thurberi × G. hirsutum) and Gossypium barbadense have been widely used in upland cotton breeding [25, 30]. First, given the knowledge that G. thurberi is a diploid D-genome cotton, we can safely assume that any introgressions from it will be vastly more likely into the upland cotton Dt subgenome (i.e., than the At subgenome). Given that the introgression-of-interest in our study located on chromosome A10 of the At subgenome, we made the assumption that G. thurberi was exceedingly unlikely to be the donor for A10 introgressions, thus narrowing the field of donors to G. arboreum and G. barbadense. Following a recent reported method for identification introgression [26, 31], we separately analyzed the introgressed fragments from G. arboreum and from G. barbadense. Fortunately, this analysis indicated that SHZDI1 shared overlap with a fragment from G. aboreum and shared no overlap with G. barbadense (Additional file 2: Table S8), findings supporting that the introgressions are from G. arboreum. A total of 36 accessions carried introgressions were identified in the GWAS panel. To further confirm the above results, we identified the SNPs between a core collection of G. arboreum germplasm resources generated by our previous report [32] with the GWAS panel on A10 (i.e., SNP-A10). The genotypes in the SHZDI1/SHZDP2/AYDP1 locus can be divided into 2 haplotypes (Hap. 1 and Hap. 2) in G. arboreum and into 5 haplotypes (Hap. A to Hap. E) in G. hirsutum (Fig. 1c). We found that Hap. B is highly similar to Hap.1, indicating Hap. B introgressed from G. arboreum. We next constructed a SNP-based phylogenetic tree using the SNPs from the putatively introgressed genomic region (highlighted in the blue box in Fig. 1c). The phylogenetic tree indicated that accessions carrying Hap. B were clustered with accessions carrying Hap. 1, further supporting that Hap. B is introgressed from G. arboreum (Fig. 1e and Additional file 1: Fig. S2a). In previous studies, G. arboreum possesses Verticillium dahliae, Fusarium oxysporum vasinfectum, and cotton leaf curl virus resistance characteristics unavailable in the tetraploid cultivated cotton gene pool [33]. Compared with accessions carrying Hap. 1, the accessions carrying Hap. 2 demonstrated no difference in VW resistance (Fig. 1d). However, when compared with G. hirsutum, accessions carrying either Hap. 1 or Hap. 2 were more tolerant to V. dahliae, indicating that G. arboreum is more tolerant to V. dahliae than G. hirsutum (Fig. 1d). We also observed that the accessions carrying Hap. B were more tolerant than those carrying Hap. A, C, D, or E (Fig. 1d and Additional file 1: Fig. S1f). Our results indicated that the introgression from G. arboreum can increase the tolerance of G. hirsutum to V. dahliae. Given that only 8.5% of accessions in the core collection of G. hirsutum germplasm carry these introgressed fragments, this favorable allele has not been widely used in cotton breeding and could be an important resource for improving VW resistance in cotton.

A large TIR-NBS-LRR gene cluster on SHZDI1 contributes to Verticillium wilt resistance

The genomic spectrum of SHZDI1/SHZDP2/AYDP1 ranges from 112.7 to 113.75 Mb on A10 (Fig. 2a). To find the corresponding homologous section on G. arboreum, we extracted 2 Mb genomic sequences (A10:112–114 Mb) of TM-1 and cut them into 1 kb fragments to blast with G. arboreum genomic sequence. The SHZDI1 locus in G. hirsutum corresponded to 132.45 to 133.60 Mb region on Chr10 of G. arboreum (Fig. 2a, b). A total of 74 genes locate in this region, of which 42 genes belong to the Toll/interleukin 1 receptor (TIR) NLR (TNL) family or its truncated genes that play an integral role in the host immune system by triggering ETI response [2, 34] (Additional file 2: Table S4 and S5).

Fig. 2

Identification of candidate genes associated with Verticillium wilt resistance on introgressed fragments (SHZDI1/SHZDP2/AYDP1). a Manhattan plot for AYDP1 on chromosome A10 and haplotype blocks (A10:112.511–113.995 Mb) around SHZDI1/SHZDP2/AYDP1 (A10:112.70–113.75 Mb) was estimated using pairwise LD correlations (R2). b Collinearity analysis of A10 of G. hirsutum and Chr10 of G. arboreum is shown in mauve. SHZDI1/SHZDP2/AYDP1 (A10:112.70–113.75 Mb) in G. hirsutum corresponded to 132.45 to 133.60 Mb on Chr10 of G. arboreum. Corresponding homology to SHZDI1/SHZDP2/AYDP1 is marked on the purple line. c Expression patterns of 4 of 14 candidate genes at 0, 3, 9, 24, 48, and 72 h post-inoculation (hpi) with V. dahliae strain Vd080 in introgressed Hap. B (ZZ2), measured by qRT-PCR. d–f Silencing of candidate genes in cotton seedlings inoculated with Vd080. Plant wilt phenotype and stem browning symptoms were photographed at 3 weeks post-inoculation. Percentages were counted for each disease grade (increasing disease severity from grade 0 to 4). DI was the mean value of three independent experiments and was measured at 21 days post-inoculation (dpi). TRV::00 was used as a negative control, and one-way ANOVA was used in the statistical analysis, and means labeled with different letters indicate significant difference at α = 0.05. Error bars represent the SD of three biological replicates

Since inducible defense signaling pathways play important roles in cotton defense against V. dahliae [10, 35], we set to identify candidate genes within the introgressed fragments (Additional file 1: Fig. S4a and Additional file 2: Table S6). By using the mRNA-seq data from plants under V. dahliae treatment, we identified 13 genes (CG1 to CG13) that were categorized as differentially expressed genes (DEGs) (adjusted P < 0.01 and absolute fold-change > 2). The expression levels of these genes were further confirmed by quantitative PCR (Fig. 2c and Additional file 1: Fig. S4e), which showed that CG01 had the highest upregulation with a 4.92-fold-change (Fig. 2c).

N-terminal region of functional TNLs, which contains a TIR domain, executes cell death response upon effector recognition, of which truncated TIR domain alone can signal effector-independent cell death response when expressed ectopically in tobacco leaves [36,37,38]. We annotated TIR domain in the 42 truncated or typical TNLs, and 28 genes encoding TIR domain were identified (Additional file 1: Fig. S4b and Additional file 2: Table S5). Next, the TIR domain from 28 genes as well as the positive (RPP1_NdATIR) and negative controls (empty vector or buffer) were transiently expressed in tobacco leaves respectively. It is interesting that only the truncated TIR domain from GA10G321300 (CG01), GA10G322400 (named CG14 but not belong to 13 DEGs), and the positive control triggered effector-independent cell death, which suggests that CG01 and CG14 have potential functions in ETI signaling pathways (Additional file 1: Fig. S4b).

To validate the potential function of the 14 candidate genes (CG01-CG14), we performed virus-induced gene silencing (VIGS) experiments. qPCR results showed that CG14 had no detectable expression, and the expression of the other 13 candidate genes was reduced significantly in cotton (Additional file 1: Fig. S4c). Plants carrying TRV::CG01 or TRV::CG04 displayed typical symptoms of VW (Fig. 2d-f). Considering the possibility of off-target silencing, the expression of two potential off-target sites with the highest homology were measured by qPCR; the transcripts of these genes were not affected in VIGS experiments (Additional file 1: Fig. S4d). Additionally, the silencing of CG01 made the introgressed resistant accession B061 (ZZ2) susceptible to V. dahliae with a highest disease index of 42.1 (Fig. 2f).

In whole, functional validation results indicate that the VW resistance mediated by introgressed fragments was associated with multiple genes in this region. Given the highly upregulated expression after inoculation, potential functions in ETI response, and the severe knockdown symptoms, CG01 was identified as one of the major VW-associated genes in this candidate region. Hereafter, we will refer to CG01 as Resistant to Verticillium dahliae 1 in Gossypium hirsutum (GhRVD1).

Allelic variation in TIR domain attenuates GhRVD1-mediated VW resistance

GhRVD1 was cloned from Hap. B (B014, ZZ2, and F045) and Hap. E (B009, B010, and ZM24) for sequence analysis. As expected, the sequences of GhRVD1 were highly divergent between Hap. B and Hap. E, with 22 nonsynonymous mutations out of 26 SNPs that were named GhRVD1_R (resistant) and GhRVD1_S (susceptible), respectively (Fig. 3a). Three domains, including the TIR domain (PFAM01582), NB-ARC domain (CL26397), and LRR domain (PFAM12799), were identified in GhRVD1 (Fig. 3a and Additional file 1: Fig. S5a). The two haplotypes of GhRVD1 were overexpressed in Arabidopsis thaliana to generate the GhRVD1_R and GhRVD1_S overexpression lines (Additional file 1: Fig. S5c). GhRVD1_R overexpression lines demonstrated enhanced resistance to V. dahliae, compared with WT and GhRVD1_S lines. Interestingly, GhRVD1_R lines demonstrated an obvious dwarfing phenotype (Additional file 1: Fig. S5c). Subsequently, two artificial transcripts of a combination of GhRVD1_S and GhRVD1_R, GhRVD1_R/S and GhRVD1_S/R, were also overexpressed in A. thaliana. GhRVD1_R/S was developed by replacing the TIR1-TIR2 (1-354AA) from GhRVD_S with that of GhRVD_R. Likewise, GhRVD1_S/R was developed by replacing the TIR1-TIR2 from GhRVD1_R with that of GhRVD_S. We found that overexpression of GhRVD1_R/S, like overexpression of GhRVD1_R, conferred A. thaliana VW resistance, while overexpression of GhRVD1_S/R showed susceptible phenotyping when compared with transgenic A. thaliana carrying GhRVD1_R (Additional file 1: Fig. S5e-g). The results indicated that the first 10 nucleotide substitutions at the position of TIR domain in the background of GhRVD1_R attenuated the acquired V. dahliae resistance.

Fig. 3

Highly divergent GhRVD1 between Hap. B and Hap. E. a An illustration of a 26 base point mutation in the CDS, and 22 nonsynonymous mutations in amino acid sequences among GhRVD1_R and GhRVD1_S; R and S are shown in green and black words, respectively. Mutation positions are marked on the top and bottom. b Phenotype of GhRVD1_R and GhRVD1_S overexpression cotton lines. c Transcript detection of both GhRVD1 genotypes using semi-quantitative PCR with universal primers; GhHiston3 was used as an internal control. d Phenotypes with a significant reduction of plant height and internode distances. Different letters indicate significant difference at α = 0.05 level via one-way ANOVA analysis. e Trypan blue staining to detect dead cells (top), and DAB staining to detect accumulation of ROS species (bottom) in 10-week-old sample leaves. f SA content of 10-week-old sample leaves via HPLC. g Expression of SAR markers (e.g., NPR1 and PR1) and HR markers (e.g., HIN1 and HSR203J) in GhRVD1 overexpression lines from leaves collected from 10-week-old plants. h Plant wilt phenotype, leaves wilt phenotype, and stem browning symptoms were photographed at 3 weeks post-inoculation. i DI was the mean value of three independent experiments and was measured at 21 dpi

Overexpression of GhRVD1 shows constitutive defense activation in cotton

We also overexpressed the two GhRVD1 haplotypes in cotton and found that 35S::GhRVD1_R induced an obvious dwarfing phenotype with a significant reduction of plant height by 61.3% and of internode distances by 40.2% (Fig. 3b–d). Dwarfing mutants have been found in cotton as well as many other plants, and brassinolide as well as gibberellin (GA) are the two major hormones associated with dwarfing [39,40,41]. Moreover, hyperactivity of plant innate immune receptors often cause ectopic defense activation, also called autoimmunity, which can manifest as severe growth retardation and spontaneous lesion formations. This phenotype is most widely studied in A. thaliana snc1 dwarfing mutant, a gain-of-function NLR mutant with constitutive defense activation [42,43,44]. Similar to the autoimmunity, overexpression of GhRVD1 induced a constitutive accumulation of reactive oxygen species (ROS) and salicylic acid (SA) and caused spontaneous lesions in 35S::GhRVD1_R (Fig. 3e, f). Hypersensitive responses (HR) induced by R genes are accompanied by several active physiological responses that restrict pathogen colonization and activate the expression of HR markers such as HIN1 and HSR203J [43, 45], both of which were significantly upregulated in all GhRVD1 overexpression lines. Moreover, the expression of HSR203J and HIN1 in 35S::GhRVD1_R was significantly enhanced compared to 35S::GhRVD1_S (Fig. 3g). The R gene-mediated HR also often triggers a secondary resistance response known as systemic acquired resistance (SAR), which is characterized by upregulation of NPR1 and PR1 [43], and we found that the expression of NPR1 and PR1 in GhRVD1 overexpression lines was significantly enhanced compared to WT (Fig. 3g). These results support that the accumulation of GhRVD1 transcripts, especially of the GhRVD1_R haplotypes, result in HR responses and an active SAR signaling, thus conferring the observed autoimmune phenotypes.

Subsequently, 35S::GhRVD1_S, 35S::GhRVD1_R, and ZM24 (negative control) plants were infected with Vd080 to evaluate the gene function. Different from the case for A. thaliana, overexpression GhRVD1_R or GhRVD1_S can both enhance resistance to VW in G. hirsutum, but a more significant enhancement of VW resistance was observed in GhRVD1_R overexpression lines (Fig. 3h, i).

Allelic variation in the TIR domain affects autoactivity of GhRVD1

Plant NLRs characterized by a multi-domain architecture consisting of either an N-terminal coiled-coil (CC) or TIR domain, a central NBS domain, and C-terminal LRR domain [46]. Typical TNLs usually have one TIR domain in the N-terminal of a protein. However, unlike the typical TNLs, GhRVD1 has two consecutive TIR domains (PF01582) located in the first 354 amino acids that share 56.4% amino acid similarity with each other, which arose by gene amplification inside cotton (Fig. 3a and Additional file 1: Fig. S6a). We screened the identified TIR domains for a more detailed comparison [36]. The first one located in the first 169 amino acids of the N-terminal, and the second TIR domain located adjacent to the first one (from 170 to 354 AA). These two TIRs from GhRVD1_R were designated as R-TIR1 and R-TIR2. Likewise, the dual-TIR domains from GhRVD1_S were named S-TIR1 and S-TIR2 in the same pattern (Fig. 3a and Additional file 1: Fig. S5b). In whole, this 354 AA region includes the TIR1 and TIR2 domains, designated as “TIR1-TIR2_R” or “TIR1-TIR2_S”.

TIR domains with known 3D structures consist of a five-stranded parallel β-sheet (βA–βE) surrounded by five α-helical regions [36]. Protein modeling with Phypre2 at 90% accuracy indicated that TIR1-TIR2_R and TIR1-TIR2_S contain most of the secondary structures in known 3D structures, except for the βB-sheet of TIR1-TIR2 and βE-sheet for TIR1 (Additional file 1: Fig. S5b). Moreover, two haplotypes of TIR1-TIR2 were detected with 10 nonsynonymous mutations in amino acid sequences (Fig. 3a). The combined mutations are predicted to result in the disappearance of αE1-helices and completely change the 3D structure of TIR1-TIR2 (Additional file 1: Fig. S5b and 6b).

In previous experiments, the autoactivity of TIR1-TIR2_R (CG01TIR) was first identified in cotton NLRs (Additional file 1: Fig. S4b), and we further investigated the divergence of autoactivity between R and S haplotypes using fusion proteins and a series of mutations. We found the intensity of cell death caused by TIR1-TIR2_S (SS) is significantly weaker than that caused by TIR1-TIR2_R (RR). To confirm whether a single TIR or both TIRs were responsible for the difference in cell death intensity, we artificially created two mutants: TIR1-S-TIR2-R (SR) and TIR1-R-TIR2-S (RS) in TIR1-TIR2_R background. When TIR1-R was mutated to TIR1-S by mutating 5 nucleotides, the cell death intensity of the first mutant SR has no difference with that of RR. However, when TIR2-R was mutated to TIR2-S by mutating 5 nucleotides, the cell death intensity of the second mutant RS significantly decreased (Fig. 4a,b and Additional file 1: Fig. S7a). Trypan blue detection of death cells also showed that RR and SR induced higher intensity of cell death compared with both SS and RS (Additional file 1: Fig. S8). The firefly Luciferase (LUC) reporter gene system has been a powerful tool to monitor the cell viability through quantifying the fluorescence activity in many studies [47,48,49,50]. Here, we also co-expressed the candidate proteins and LUC to perform the cell death assay. Both the fluorescence signal and luciferase activity again confirm that RR and SR induced stronger intensity of cell death than both SS and RS (Fig. 4d,e and Additional file 1: Fig. S6h). Viewed together, our results suggest that the HR response’s intensity enhanced in a manner depended on the mutations between the TIR2-R and TIR2-S.

Fig. 4

Nonsynonymous SNP sites of 319 and 326 in TIR1-TIR2 caused differences in HR strength between the two haplotypes, and the minimum functional region is TIR1-TIR2. a,b HR phenotypes associated with chimeras or site-directed mutants of TIR1-TIR2_R and TIR1-TIR2_S alleles. RS and SR are artificial constructs of a combination of R-TIR1-S-TIR2 and S-TIR1-R-TIR2, respectively. Site-directed mutants were converted to S alleles at the nonsynonymous SNP sites in a TIR1-TIR2_R background. c Evaluation of protein expression using a Flag antibody, and immunoblotting of plant actin with α-actin was used as a loading control. d,e LUC activity within the region co-transformed with a LUC plasmid and TIR1-TIR2 variants. Fluorescence intensity was captured and data represent means ± SE of three independent experiments. Different letters indicate significant difference at α = 0.05 level via one-way ANOVA analysis. f,g Cell death phenotype in cotton protoplast for cell viability assays. Protoplast preparation and construct transformation used 14-day-old plant leaves. Fluorescence intensity was captured. Three independent experiments are represented. h Schematic diagram of the GhRVD1 domain structure. Individual domains are presented in a colored box, and the boundaries of truncated sequences are marked on the top. The numbers and names of truncated sequences are marked on the right and left, respectively. i HR phenotypes of truncated derivatives of GhRVD1_R. j Evaluation of protein expression using a Flag antibody, immunoblotting of plant actin with α-actin was used as a loading control. k LUC activity within the region co-transformed with a LUC plasmid and a series of truncated derivatives of GhRVD1. Numbers on the top represent corresponding truncated sequences

Following this reasoning, K192N, K299Q, K319T, K326E, K319T/K326E, and L344F were mutated at the nonsynonymous SNP sites of TIR2 in a TIR1-TIR2_R background (Fig. 4a) and were expressed ectopically in tobacco leaves. The results demonstrate that amino acid substitutions at positions 319 or 326 attenuated the HR response (Fig. 4b and Additional file 1: Fig. 7a), which locate in α-helix E1 and α-helix E2 of TIR2, respectively (Additional file 1: Fig. S6b), but the strength of cell death was not further attenuated in K319T/K326E. Mutations at these two amino acids caused differences in HR strength between the two haplotypes of TIR1-TIR2. Protein immunoblotting demonstrated that all fusion or mutant proteins were properly expressed in tobacco leaves (Fig. 4c). Subsequently, trypan blue assay (Additional file 1: Fig. S8) and fluorescence intensity (Fig. 4d,e) again supported that amino acid substitutions at position 319 or 326 of TIR2, respectively, attenuated the HR response.

To verify that TIRs can result in cell death in cotton, TIR1-TIR2 variants and 35S::Luciferase were simultaneously expressed in cotton protoplasts. The decrease of the fluorescence signal indicates that TIR1-TIR2 also caused cell death in cotton protoplasts, and the cell death caused by SS and K319T/K326E was significantly weaker than that of RR (P < 0.01) (Fig. 4f,g). Our results indicate that the two haplotypes of TIR1-TIR2 from GhRVD1 can induce cell death, and nonsynonymous SNP sites 319 and 326 of TIR1-TIR2 play an important role in controlling the intensity of cell death.

Defining the functional boundaries of the GhRVD1 TIR domain

To clarify the functional boundaries of the TIR domain, a series of truncated GhRVD1 derivatives were constructed and transiently expressed in tobacco leaves (Fig. 4h). We found that NBS-ARC sequences connected after the dual-TIR domains can significantly inhibit cell death, which was similar to the previous report assessing RPP1_NdA [37]. Notably, neither TIR1 nor TIR2 alone could induce cell death, and the smallest fragment that induced cell death was amino acid 1–354 of GhRVD1, which precisely covered the full length of the TIR1 and TIR2 domain sequences (Fig. 4h–j and Additional file 1: Fig. S7b). With the exception of LRR-L, all proteins were detectable by protein immunoblotting, although the steady-state accumulation of each protein differed substantially (Fig. 4j). Our results suggest that this 354 AA region, comprising “TIR1-TIR2,” is the minimal sequence required for autoactivity and was inhibited by sequence ligation of the NB-ARC domain.

Mutations at the predicted DE interface disrupt autoactivity and self-association of TIR1-TIR2

Plant TNL activation leads to homodimerization of the intracellular TIR domain and initiates the downstream signaling pathway [2, 51]; this self-association is also present in the effector-independent cell death of the TIR domain autoactivation phenotype [1, 36, 46]. In a previous study, predicted αA- and αE-helices (AE) and DE interfaces were found to be involved in the homodimerization of the plant TIR domain [1, 36, 52]. Based on homologous alignment with the well-characterized TIR domains of L6 [36] and RPS4 [