Jurassic NLR: Conserved and dynamic evolutionary features of the atypically ancient immune receptor ZAR1

Abstract

Plant nucleotide-binding leucine-rich repeat (NLR) immune receptors generally exhibit hallmarks of rapid evolution, even at the intraspecific level. We used iterative sequence similarity searches coupled with phylogenetic analyses to reconstruct the evolutionary history of HOPZ-ACTIVATED RESISTANCE1 (ZAR1), an atypically conserved NLR that traces its origin to early flowering plant lineages ∼220 to 150 million yrs ago (Jurassic period). We discovered 120 ZAR1 orthologs in 88 species, including the monocot Colocasia esculenta, the magnoliid Cinnamomum micranthum, and most eudicots, notably the Ranunculales species Aquilegia coerulea, which is outside the core eudicots. Ortholog sequence analyses revealed highly conserved features of ZAR1, including regions for pathogen effector recognition and cell death activation. We functionally reconstructed the cell death activity of ZAR1 and its partner receptor-like cytoplasmic kinase (RLCK) from distantly related plant species, experimentally validating the hypothesis that ZAR1 evolved to partner with RLCKs early in its evolution. In addition, ZAR1 acquired novel molecular features. In cassava (Manihot esculenta) and cotton (Gossypium spp.), ZAR1 carries a C-terminal thioredoxin-like domain, and in several taxa, ZAR1 duplicated into 2 paralog families, which underwent distinct evolutionary paths. ZAR1 stands out among angiosperm NLR genes for having experienced relatively limited duplication and expansion throughout its deep evolutionary history. Nonetheless, ZAR1 also gave rise to noncanonical NLRs with integrated domains and degenerated molecular features.

Figure 1.

Comparative sequence analyses identify and classify ZAR1 sequences from angiosperms. A) Workflow for computational analyses in searching for ZAR1 orthologs. We performed TBLASTN/BLASTP searches and subsequent phylogenetic analyses to identify ZAR1 ortholog genes from plant genome/proteome datasets. B) ZAR1 forms a clade with 2 closely related sister subclades. The phylogenetic tree was generated in MEGA7 by the neighbor-joining method using NB-ARC domain sequences of ZAR1-like proteins identified from the prior BLAST searches and 1,019 NLRs identified from 6 representative plant species, taro, stout camphor, columbine, tomato, sugar beet, and Arabidopsis. Red arrowheads indicate bootstrap support > 0.7 and is shown for the relevant nodes. The scale bar indicates the evolutionary distance in amino acid substitution per site.

Figure 2.

The ZAR1 gene is distributed across angiosperms. The phylogenetic tree was generated in MEGA7 by the neighbor-joining method using full-length amino acid sequences of 120 ZAR1 orthologs identified in Fig. 1. Red triangles indicate bootstrap support > 0.7. The scale bar indicates the evolutionary distance in amino acid substitution per site.

Figure 3.

ZAR1 orthologs carry conserved sequence patterns required for Arabidopsis ZAR1 resistosome function. A) Schematic representation of the Arabidopsis ZAR1 protein highlighting the position of conserved sequence patterns across ZAR1 orthologs. Consensus sequence patterns were identified by MEME using 117 ZAR1 ortholog sequences. Raw MEME motifs are listed in Supplemental Table S1. Red asterisks indicate residues functionally validated in Arabidopsis ZAR1 for NBD–NBD and ZAR1–RLCK interfaces. B) Conservation and variation of each amino acid among ZAR1 orthologs across angiosperms. Amino acid alignment of 117 ZAR1 orthologs was used for conservation score calculation via the ConSurf server (https://consurf.tau.ac.il). The conservation scores are mapped onto each amino acid position in Arabidopsis ZAR1 (NP_190664.1). C, D) Distribution of the ConSurf conservation score on the Arabidopsis ZAR1 structure. The inactive ZAR1 monomer is illustrated in cartoon representation with domain architecture C) and conservation score D). Major 5 variable surfaces (VS1 to VS5) on the inactive ZAR1 monomer structure are described in gray dot or black boxes in panel B or D, respectively.

Figure 4.

ZAR1 orthologs across angiosperms display multiple conserved surfaces on the resistosome structure. Distribution of the ConSurf conservation score was visualized on the inactive monomer A), active monomer B), and resistosome C) structures of Arabidopsis ZAR1. Each structure and cartoon representation are illustrated based on the conservation score shown in Fig. 3.

Figure 5.

ZRK gene clusters occur in A. coerulea and C. micranthum. A) The phylogenetic tree was generated in MEGA7 by the neighbor-joining method using full-length amino acid sequences of 39 ZRK proteins. Red triangles indicate bootstrap support > 0.7. The scale bar indicates the evolutionary distance in amino acid substitution per site. B) Schematic representation of the ZRK gene clusters on an A. coerulea (columbine) contig and a C. micranthum (Stout camphor) scaffold.

Figure 6.

ZRK family proteins positively regulate A. coerulea AcZAR1 and C. micranthum CmZAR1 autoimmune cell death in N. benthamiana. A, C) Cell death observed in N. benthamiana after expression of ZAR1 mutants with or without wild-type ZRKs. N. benthamiana leaf panels expressing wild-type NbZAR1 (NbZAR1^WT), NbZAR1^D481V (ZAR1^D481V), AcZAR1^D489V (AcZAR1^DV), and CmZAR1^D488V (CmZAR1^DV) with or without wild-type ZRKs, were photographed at 4 d after agroinfiltration. B, D) Violin plots show AcZAR1 and CmZAR1 cell death intensity scored as an HR index based on 12 and 9 replicates (different leaves from independent plants) in 2 independent experiments. Statistical differences among the samples were analyzed with Tukey's HSD test (P < 0.01).

Figure 7.

Cassava and cotton ZAR1-ID carry an additional Trx domain at the C terminus. A) Schematic representation of NLR domain architecture with C-terminal Trx domain. B) Description of Trx domain sequences on amino acid sequence alignment. Cassava XP_021604862.1 (MeZAR1) and cotton KAB1998109.1 (GbZAR1) were used for MAFFT version 7 alignment as representative ZAR1-ID. Arabidopsis ZAR1 (AtZAR1) was used as a control of ZAR1 without ID.

Figure 8.

ZAR1-SUB has emerged early in eudicots and diverged at MADA motif sequence. The phylogenetic tree was generated in MEGA7 by the neighbor-joining method using full-length amino acid sequences of 120 ZAR1, 129 ZAR1-SUB, and 11 ZAR1-CIN identified in Fig. 1. Red triangles indicate bootstrap support > 0.7. The scale bar indicates the evolutionary distance in amino acid substitution per site. NLR domain architectures are illustrated outside of the leaf labels: MADA is red, CC is pink, NB-ARC is yellow, LRR is blue, and other domain is orange. Black asterisks on domain schemes describe truncated NLRs or potentially misannotated NLR.

Figure 9.

Conserved sequence distributions in ZAR1-SUB and ZAR1-CIN. A) Schematic representation of the ZAR1-SUB protein highlighting the position of the representative conserved sequence patterns across ZAR1-SUB. Representative consensus sequence patterns identified by MEME are described on the scheme. Raw MEME motifs are listed in Supplemental Tables S2 and S3. B) Conservation and variation of each amino acid among ZAR1-SUB and ZAR1-CIN. Amino acid alignment of 129 ZAR1-SUB or 8 ZAR1-CIN was used for conservation score calculation via the ConSurf server (https://consurf.tau.ac.il). The conservation scores are mapped onto each amino acid position in queries XP_004243429.1 (ZAR1-SUB) and RWR85656.1 (ZAR1-CIN), respectively. C) Schematic representation of the ZAR1-CIN protein highlighting the position of the representative conserved sequence patterns across 8 ZAR1-CIN. Raw MEME motifs are listed in Supplemental Tables S4 and S5.

Figure 10.

Coevolution of ZAR1 and ZRK genes in angiosperms. A) We propose that the ancestral ZAR1 gene has emerged ∼220 to 150 million yrs ago (Mya) before monocot and eudicot lineages split. ZAR1 is a widely conserved CC-NLR in angiosperms, but it is likely that ZAR1 was lost in the monocot lineage, Commelinales. A sister clade paralog ZAR1-SUB has emerged early in the eudicot lineages and may have been lost in Caryophyllales. Another sister clade paralog ZAR1-CIN was duplicated from the ZAR1 gene and expanded in the Magnoliidae C. micranthum. Trx domain integration to C terminus of ZAR1 has independently occurred in few rosid lineages. B)ZAR1 has coevolved with partner ZRK gene for pathogen effector recognition since the Jurassic era. During the coevolution, ZRKs diversified to catch up with fast-evolving effectors.