Phylogenomic analysis reveals exceptions to the co-evolution of ZAR1 and ZRK immune gene families in plants

Abstract

Background

HOPZ-ACTIVATED RESISTANCE 1 (ZAR1) is a nucleotide-binding leucine-rich repeat (NLR) protein functioning as a recognition hub to initiate effector-triggered immunity against bacterial pathogens. To initiate defense, ZAR1 associates with different HOPZ-ETI-DEFICIENT 1 (ZED1)-Related Kinases (ZRKs) to form resistosomes to indirectly perceive effector-induced perturbations. Few studies have focused on the phylogenomic characteristics of ZAR1 and ZRK immune gene families and their evolutionary relationships. To trace the origin and divergence of ZAR1 and ZRK immune gene families across the plant kingdom, we performed phylogenomic analyses using an extended set of plant genomes.

Results

Genome-wide identification of ZAR1 and ZRK immune gene families by blast similarity searches combined with phylogenetic analysis showed that these two gene families have experienced frequent gene losses in massive lineages. Gene distribution patterns across the plant kingdom revealed that ZAR1 and ZRK emerged after the divergence of most angiosperms from Amborella and before the split of magnoliids, monocots, and eudicots. Co-occurrence of ZAR1-A and ZRKs was found in various plant species belonging to different angiosperm orders, but both genes were found to be absent in chlorophyta, bryophytes, lycophytes, ferns, and gymnosperms. We also detected a large number of concerted gene losses in angiosperms, especially within the orders Fabales, Cucurbitales, Asterales, and Apiales. All analysed monocot genomes thus far examined, except for the aroid Colocasia esculenta, were previously reported to lack both ZAR1-A and ZRKs. Here we now report other exceptions on the concerted ZAR1-A–ZRKs presence-absence pattern within several early diverging monocot lineages, including the genome of Acorus tatarinowii—a species representing the first branching monocot lineage. We also revealed strong variation in ZAR1-A–ZRKs co-occurrence within the asterid order Ericales, suggesting patterns of de-coevolution in angiosperms. Our research further shows that both gene families experienced significant diversification through various duplication events. Additionally, their evolutionary paths have been shaped by frequent gene losses and lineage-specific transposition.

Conclusion

This study provides novel findings on the evolution of ZAR1 and ZRK immune gene families across a wide range of plant species, suggesting that more potential exceptions can be expected when expanding the list of sequenced genomes from distinct orders. Our results provide new hypotheses about the origin and diversification of these critical immune genes for future functional studies.

Background

Phenotypic diversity is prevalent in various angiosperm lineages, and each lineage has its own trait characteristics [1]. This phenomenon is mainly caused by patterns of gene diversification. It underpins the ability of species to evolve to resist diverse environmental stresses and is thought to be largely driven by frequent gene duplication and consequent functional divergence [2]. Polyploidization events, including whole-genome duplication (WGD) and whole-genome triplication (WGT), are considered to have a prominent effect on gene copy number [3]. Ancestral polyploidization is often followed by diploidization events, which largely reorganize genome structure by repetitive DNA loss, chromosomal rearrangements, and multiple gene losses [4]. Besides, gene duplication can also occur on a local scale through multiple processes, i.e. tandem, proximal, and dispersed duplication [5]. Gene duplicates derived from both genome-wide and local duplication processes can be preserved by neo- and sub-functionalization, leading to gene diversification that can influence trait evolution [6]. Preservation of genomic context during evolution is known as synteny, which is regarded as a sign of functional conservation of genes [7]. Visualization of genomic rearrangements and gene synteny can be used to decipher the dynamics of gene family evolution.

Plant effector-triggered immunity (ETI) is mediated by nucleotide-binding leucine-rich repeat (NLR) proteins, which recognize the activity of effector proteins secreted by invasive pathogens and elicit a robust plant defense, often accompanied with a hypersensitive response (HR) [8]. Recognition of effectors by NLR proteins is generally indirect due to the co-evolution of host–pathogen interactions [9]. This recognition in some cases is achieved through a mechanism known as the guard model, in which pathogen effectors modify host target proteins (guardees) to activate NLR-mediated ETI. Alternatively, NLR proteins can perceive effector-triggered modifications of host decoy proteins via structures similar to true effector targets, thus trapping pathogen effectors into plant recognition system and inducing downstream immune responses [9,10,11]. Indirect effector recognition allows plants to obtain a broad spectrum of pathogen recognition with limited NLR proteins [9]. A canonical example is HOPZ-ACTIVATED RESISTANCE 1 (ZAR1), a CC-type NLR that has been regarded as a recognition hub capable of sensing multiple effectors of different Arabidopsis pathogens, such as the effectors HopZ1a, HopF2a, HopBA1, HopX1, and HopO1 from Pseudomonas syringae [12,13,14], and the effector AvrAC from Xanthomonas campestris [15]. This striking diversity in the recognition of effectors by ZAR1 occurs through interactions with members of the receptor-like cytoplasmic kinase XII-2 (RLCK XII-2) subfamily, also known as HOPZ-ETI-DEFICIENT 1 (ZED1)-associated kinases (ZRKs) [12,13,14,15].

The Arabidopsis RLCK XII-2 subfamily consists of 13 ZRKs, several of which are known to play roles in ZAR1-mediated immunity towards bacterial pathogens. Recognition of the P. syringae effector HopZ1a in Arabidopsis was found to rely on a preformed complex consisting of ZAR1 and ZED1 (also known as ZRK1) [14]. The acetylation of AtZED1 by HopZ1a is shown to activate AtZAR1, thus triggering downstream immune responses [14]. Recognition of another P. syringae effector, HopF2a/HopF1r, requires the involvement of ZRK3 to induce ZAR1-mediated immunity. More recently, ZRK1, ZRK2, and ZRK3 were found to interact with ZAR1 to regulate the recognition of P. syringae effectors HopX1, HopBA1, and HopO1, respectively—further demonstrating the central role of ZAR1 in immune activation to ward off infection by Pseudomonas pathogens [16]. In addition, a similar recognition mechanism has been assigned to the recognition of the X. campestris effector AvrAC [15]. The ZAR1-ZRK1 complex recruits the decoy protein PBL2 to participate in AvrAC recognition. Uridylylation of this decoy protein by AvrAC promotes the binding of PBL2 to ZRK1 in the ZAR1-ZRK1 complex, thereby resulting in the formation of a wheel-like pentamer [17, 18]. This so-called pentameric resistosome was suggested to induce cell death by disrupting membrane integrity [17, 19]. Beside these aforementioned functional studies in the model plant Arabidopsis, it was shown that ZAR1 of Nicotiana benthamiana (NbZAR1) interacts with the ZRK family member JIM2 to perceive the effector XopJ4 from Xanthomonas perforans [20], showing that ZAR1-ZRK interactions are not Brassicaceae-specific.

A multitude of studies have pinpointed the conserved architecture of ZAR1 in a number of dicotyledonous plants [20,21,22,23]. Phylogenomic analysis of ZAR1 across 88 angiosperm species revealed that ZAR1 first appeared ~ 220 to 150 million years ago (Mya) [24]. A recent study that made use of 72 plant genomes proposed that the ZAR1-ZRK resistosome emerged in early angiosperms and experienced concerted gene losses across multiple lineages during the co-evolution of ZAR1 and ZRKs [25].

Phylogenomic studies with incomplete taxon sampling and a limited number of representative genomes cannot often provide a global picture of gene family evolution. To gain further knowledge on how ZAR1 and ZRK immune gene families evolved, we incorporated more representative species in each category of plants to perform a more comprehensive phylogenetic analysis. Our results confirmed the early origin and co-evolution pattern of these two immune gene families as suggested in previous studies. It corroborates that frequent gene loss events have occurred in multiple plant lineages. Our analyses, however, revealed extensive variation in the presence-absence patterns of ZAR1-A–ZRK within early-diverging monocot lineages, including various aroid species (Araceae) and Acorus tatarinowii, a species belonging to Acorales—the sister lineage to all the other extant monocot plants. This implies that other early diverging monocot taxa may also contain ZAR1-A and/or ZRKs. Clear exceptions on ZAR1-A–ZRK co-occurrence were also identified for several dicot species belonging to the order Ericales, suggesting potential divergence from the shared co-evolution mode. Moreover, our synteny analysis revealed that, next to WGD and tandem duplication, genes of both families diversified by lineage-specific transposition events. This study extends our understanding of the evolution of ZAR1 and ZRK immune gene families across the plant kingdom and provides a theoretical basis for future functional studies on plant disease resistance.

Results

Broad-scale identification of ZAR1 and ZRK genes in plants

A collection of 231 genomes from diverse plant species (Additional Table 1) [26]—ranging from green algae to eudicots—was used for extensive phylogenomic analyses of ZAR1 and ZRK immune gene families (Additional Fig. 1). A total of 254 ZAR1 and 929 ZRK homologous genes were detected in 100 and 99 plant genomes, with an average of 2.5 and 9.4 members per genome, respectively (Additional Table 2–3). Copy number variation was determined by calculating coefficients of variation (CVs)—i.e. ratios of the standard deviation (SD) to the mean, expressed as a percentage. CVs of ZAR1 and ZRK genes were estimated around 84% and 89%, respectively. These values were particularly high when compared with conserved gene families that play a role in plant development such as the CesA family, which has a CV of approximately 40% [26]. This large difference could suggest that immune genes are prone to have higher variability caused by evolutionary adaptation to biotic stresses. Remarkably, ZAR1/ZRK copy numbers were found not to be correlated with ploidy level, and only weakly correlated with the number of genome duplication events and genome sizes across the analyzed species (Additional Fig. 2–4), implying that other factors, such as tandem duplication, strongly influenced the expansion of ZAR1 and ZRK immune gene families during plant evolution.

Expansion differences in ZAR1 and ZRK gene families

We performed phylogenetic analysis to explore the evolutionary relationships of the ZAR1 and ZRK immune gene families. Identified ZAR1 genes were categorized into three closely related groups, i.e. the AtZAR1-containing clade (termed ZAR1-A), the sister clade ZAR1-B, and the first splitting branch ZAR1-C (Additional Fig. 5), which is consistent with the results of Gong et al. (2022) and Adachi et al. (2023) [24, 25]. Group ZAR1-A contains ZAR1 orthologs from magnoliids, few monocots and most eudicot clades, suggesting that the ZAR1-A group was formed before the divergence of monocots, magnoliids, and eudicots. Group ZAR1-B and ZAR1-C include ZAR1 paralogs from eudicots and magnoliids, respectively (Additional Fig. 5). Similar phylogenetic analysis provided evidence that ZRK genes emerged after the divergence of most angiosperms from the Amborella lineage (Additional Fig. 6), indicating that the origin of ZRK could be traced back to the root of monocots, magnoliids, and dicotyledons. The observation that the phylogeny of ZRKs fall into distinct lineages implies a complex evolutionary process involving extensive asymmetric expansion and loss.

To systematically investigate the expansion of ZAR1 and ZRK genes across genomes, the gene copy number belonging to each group within each species was mapped onto the corresponding leaves (tips) of the species tree (Fig. 1 and Additional Table 4). Absence of both ZAR1 and ZRK families was observed in multiple early diverging plant lineages, including chlorophyta, bryophytes, lycophytes, ferns, and gymnosperms (Fig. 1). Nevertheless, the distribution of these two immune gene families differed markedly in gene copy number after the divergence of the early branching angiosperm Amborella (Fig. 1 and Additional Fig. 7). The overall gene copy number of the ZAR1-A group was relatively low, i.e. between 1 and 4 (CV = 42.4%), while the ZAR1-B and ZAR1-C group confined to eudicots and magnoliids have relatively higher copy numbers ranging from 1–9 (CV = 74.0%) and 1–11 (CV = 98.0%), respectively. However, ZRKs were greatly expanded in angiosperms, especially in Brassicales and Fagales, with average copy numbers of 17 (CV = 59.7%) and 21 (CV = 76.8%), respectively. It should be noted that Caricaceae and Cleomaceae, two families within Brassicales, lack or have relatively low copy number of ZRKs (Fig. 1 and Additional Table 4). Exploring gene duplication patterns showed that most ZRKs in Brassicales and Fagales were derived from dispersed and/or tandem duplication (Fig. 2). Inter-species differences in expansion patterns were also obvious. For example, the allotetraploid species B. napus and its diploid progenitors B. oleracea and B. rapa have a large proportion of ZRK genes originating from WGD or segmental duplication. A total of 28 ZRK homolog pairs were detected in B. napus, of which 15 (53.6%) of these duplication events occurred between A and C sub-genomes due to the polyploid nature of B. napus (Additional Fig. 8). Similarly, the high number of ZRKs in Camelina sativa can be explained by polyploidization and WGD [27]. However, no ZRK genes generated via this pattern were present in Fagales, especially in Quercus suber that displays a comparable high copy number of ZRKs (Fig. 2).

Fig. 1

Distribution of ZAR1 and ZRK genes in 231 plant genomes. The size of each circle reflects the gene copy number per category (ZAR1-C (a), ZAR1-B (b), ZAR1-A (c), and ZRK (d)). Pink, red, and blue stars represent reported genome additions, WGD, and WGT events, respectively [28,29,30,31]

Fig. 2

Percentage of ZRK genes derived from different duplication patterns in Brassicales (A) and Fagales (B). Bar plots represent the percentage of ZRK genes derived from different duplication patterns. Copy number variation across species is indicated by the size of the turquoise-colored circles. Pink, red, and blue stars represent reported genome additions, WGD, and WGT events, respectively [28,29,30,31]. Phylogenetic relationships were extracted from Fig. 1

To trace the ancestral status of ZAR1 and ZRK immune gene families, we assessed two distinct characters for each category in the Angiosperm Phylogeny Group (APG) tree, one for the ratio of species containing ZAR1 or ZRK genes (species ratio), and the other for their total gene copy numbers. ZAR1-A and ZRK were both basically present in Malpighiales, Rosaceae, Cannabaceae, Fagales, Sapindales, Malvales, Brassicaceae, Caryophyllales, Gentianales, and Solanales, with a species ratio higher than 0.8 and species number higher than 3 (Table 1). However, the total copy numbers of ZAR1-A and ZRKs across these categories varied widely, ranging from 3 to 24 and 22 to 369, respectively. This difference suggests that these two families have different evolutionary histories. In contrast, a large number of recent losses were found within both immune gene families, especially in Fabales, Cucurbitales, Asterales, Apiales, and most monocots (Table 1). In this regard, our findings align with those of Adachi et al. [24] that were based on a much smaller set of plant genomes (88) and excluded genomes from gymnosperms, ferns, mosses, charophytes, and chlorophytes.

Table 1 Species ratio and total copy number of ZAR1 and ZRK genes within different angiosperm categories

Exceptions to the co-occurrence pattern of ZAR1-A–ZRK

To further study the co-occurrence pattern between ZAR1-A and ZRKs [25], we compared gene copy numbers of ZAR1-A and ZRKs across aforementioned plant genomes. Results showed concerted presence or absence of both genes in nearly all 231 analyzed genomes (Fig. 1 and Additional Table 4). We, however, found extensive variations in the presence-absence pattern of ZAR1-A–ZRK across several early diverging species within the monocot family Araceae (Fig. 3). The aroid Amorphophallus konjac was identified as having ZRKs but lacking ZAR1-A, whereas species within the subfamily Lemnoideae (i.e. the morphologically reduced surface-dwelling duckweeds Spirodela polyrhiza and Lemna minor) and the fully submerged marine seagrass Zostera marina were found to lack both ZAR1-A and ZRKs. This co-absence may result from significant gene loss that occurred as part of their extreme adaptation to an aquatic lifestyle [32, 33]. A consequent lack of a ZAR1-ZRK-mediated resistosome might be tolerated, given that aquatic habitats typically have fewer biotic stressors than terrestrial environments.

Fig. 3

Distribution of ZAR1 and ZRK genes in Ariods (Araceae) and related monocots. Circle sizes reflect the gene copy numbers of ZAR1-A and ZRK. Red stars represent reported WGD events [28, 30, 31]. Phylogenetic relationships were extracted from Fig. 1

Besides Colocasia esculenta (taro) [24, 25], we identified co-occurrence of ZAR1-A–ZRK in the closely related surface-dwelling species Pistia stratiotes. Unlike duckweeds, Pistia forms rosettes and has potentially retained the ZAR-ZRK resistosome due to its need for pathogen protection throughout its longer life span. More interestingly, a similar co-occurrence pattern was observed in Acorus tatarinowii, a species belonging to the first branching monocot lineage Acorales, which is sister to the Araceae clade and all other monocot lineages (Fig. 3, Additional Table 4). We also found two angiosperm species that do not fit this co-evolution pattern, i.e. Solanum melongena (eggplant) and Pyrus bretschneideri (pear) that both lack ZAR1-A but do have ZRKs (Fig. 1, Additional Table 4). Considering that other genomes from Solanales and Rosales do show a ZAR1-A–ZRK co-occurrence pattern, these two exceptions can be likely explained by errors caused by genome assembly or annotation. Exceptions that clearly do not match this co-occurrence pattern can be found when analyzing several genomes of Ericales (Additional Table 1). Here we showed that the genomes of Vitellaria paradoxa, Camellia sinensis, Actinidia chinensis, Vaccinium darrowii, and Rhododendron simsii do contain ZRKs but lack ZAR1-A, indicating potential divergence from the shared co-evolution mode between ZAR1-A and ZRKs (Fig. 4). In addition, we also investigated ZAR1-B and ZRK genes in species that lack ZAR1-A. Camellia sinensis and two aforementioned species potentially lacking ZAR1-A, i.e., eggplant and pear, do however contain two or three ZAR1 paralogs from the ZAR1-B group (Additional Table 4), indicating potential sub-functionalization. However, we did not observe any species containing both ZAR1-C and ZRK genes in the absence of ZAR1-A and ZAR1-B.

Fig. 4

Distribution of ZAR1 and ZRK genes in genomes of Ericales and related asterids. Circle sizes reflect gene copy numbers of ZAR1-C, ZAR1-B, ZAR1-A, and ZRK. Red and blue stars represent reported WGD and WGT events, respectively [28, 30, 31]. Phylogenetic relationships were extracted from Fig. 1

Proximity of ZAR1 and ZRKs in chromosome-level assembled genomes

While NLR genes usually cluster together to facilitate their rapid diversification and evolution [34, 35], the same does not apply to ZAR1 which is 183 kb from the nearest NLR in Arabidopsis [24]. Since ZAR1 is situated on the same chromosome with a ZRK cluster in Arabidopsis, and interacts with diverse ZRKs to sense different pathogenic effectors [22], we speculate that this position relationship may be advantageous.

To test this hypothesis and avoid the negative effects of low-quality genomes, we selected 40 chromosome-level assembled genomes containing both ZAR1 and ZRK genes to explore the potential genomic co-location of members of the two immune gene families (Additional Table 5). Results showed that 19 out of 40 analyzed angiosperm species (47.5%) have both ZAR1-A and/or ZAR-B genes located on the same chromosome with ZRK genes. This phenomenon can be found across diverse families, including Brassicaceae, Araceae, and Salicaceae. In total, 92.0% of ZAR1-A genes and 23.5% of ZAR1-B genes in these 19 species were found to co-situate with at least one ZRK on the same chromosome. However, most of these ZAR1-A and/or ZAR1-B genes were distantly positioned from the nearest ZRK gene (Additional Fig. 9), indicating that the proximity of ZAR1 and ZRKs occurs randomly across different genomes. This suggests that genes encoding components of protein complexes are not necessarily located in close proximity to each other on the genome. It should be noted that several Brassicaceae species with genomes assembled at chromosome-level have their ZAR1-A genes (10) located near genomic ZRK clusters, which was not observed in Salicaceae. Thus, these results do not suggest a direct relationship between gene functionality and genomic position.

Conservation and variation of gene structures

To understand how gene structures of ZAR1 and ZRKs changed over evolutionary time scales, we conducted exon–intron analysis for all retrieved genes. The gene lengths of ZAR1 homologs belonging to the three groups were much longer than what were found for ZRKs (Fig. 5A). Accordingly, CDS lengths of ZAR1 homologs were at least twice as long as those of ZRKs (Fig. 5B). In some cases, gene lengths of ZAR1-B and ZAR1-C extended to 60 kb, despite their average CDS lengths of approximately 2.5 kb, suggesting the presence of exceptionally long introns possibly resulting from transposon insertions [36]. However, no significant difference was detected in exon numbers within these two immune gene families (Fig. 5C).

Fig. 5

Gene structure variation and Ka/Ks values of ZAR1 and ZRK genes. A-C, graphical display of gene length (A), CDS length (B), and exon number (C) in all collected genomes. D, Ka/Ks ratios of ZAR1-A and ZRK in 11 representative plant genomes

We used 13 representative high-quality genomes to analyze the gene structure of ZAR1 and ZRK genes in greater detail. Since most ZAR1-C and ZAR1-B homologs were absent in these genomes, we focused primarily on the gene structure variation of ZAR1-A and ZRKs. Generally, the gene lengths of ZAR1-A varied considerably across different species, while their CDS lengths were relatively conserved (Table 2). This may be caused by insertion of introns, but more importantly due to the variable length of UTRs, since most ZAR1-A genes that contain only one exon (without introns) have a much longer gene length than CDS. In addition, ZAR1-A genes from taro (monocots), Cinnamomum kanehirae (magnoliids), and Amaranthus hypochondriacus (superasterids) have more exons than that of the other angiosperms (p < 0.05). This result may indicate that the structure of ZAR1-A genes evolved from complex to simple after its emergence. Compared to their CDS lengths, ZRKs also exhibited considerable variability in gene lengths (Table 2), emphasizing substantial differences in intron lengths among ZRKs across different plant species. We also found that the average exon number of ZRK was higher than 1 in all representative species, indicating that at least one copy of ZRK in each species had at least two exons.

Table 2 Structural features of ZAR1 and ZRK genes in 13 representative species

Purifying selection of ZAR1-A and ZRK genes

To evaluate the selection pressure underlying ZAR1-A and ZRK genes, Ka/Ks ratios of each homologous pair in 11 representative plant species were calculated using ParaAT2.0 software. The results showed that Ka/Ks ratios of all homologous pairs were found to be less than 1, but the Ka/Ks of ZAR1-A genes was significantly lower than that of ZRKs (Fig. 5D). This reveals that both ZAR1-A and ZRK genes have undergone purifying selection, but ZAR1-A experienced stronger evolutionary pressure with lower Ka/Ks in all selected species, indicating that the ZAR1-A genes are more conserved than ZRKs. Considering that an individual ZAR1-A can bind different ZRKs to form diverse functional complexes to recognize different pathogen effectors [22], we infer that this difference in conservation degree between ZAR1-A and ZRK genes is logical and understandable. This result suggests that ZRKs evolved adaptively in response to infection by various pathogens over the long evolutionary history of angiosperms.

Gene duplication and synteny network analyses

The syntenic conservation of ZAR1 and ZRK immune gene families was evaluated using phylogenetic synteny network analyses based on their genomic contexts in each genome (Additional Fig. 1). Approximately 55% of the ZAR1 (145 out of 254) and ZRK (510 out of 929) genes are contained in ZAR1 and ZRK synteny networks with 2,301 and 1,772 connections, respectively (Additional Table 6). These ratios are rather low when compared to other synteny studies on highly expanded gene families, such as APETALA2, MADS-box, and CesA/Csl [26, 28, 30]. This result indicates that the collinearity of both ZAR1 and ZRK genes is less interconnected and conserved than that of several large gene families that play roles in regulatory and developmental processes. By decomposing synteny networks into clusters of closely connected genes, 6 ZAR1 and 38 ZRK syntenic communities were identified, containing 136 (93.8%) ZAR1 and 334 (65.5%) ZRK syntelogs, respectively (Fig. 6 and Additional Fig. 10).

Fig. 6

Phylogenetic profiles of ZAR1 and ZRK communities. Rows and columns display species and syntenic communities, respectively. The number of genes (syntelogs) per species within each community is indicated by a color gradient. Pink, red, and blue stars represent reported genome additions, WGD, and WGT events, respectively [28,29,30,31]

The phylogenetic profiles of ZAR1 and ZRK syntenic communities are notably different (Fig. 6). Two larger ZAR1 communities (1 and 2) were relatively conserved, consisting of syntelogs from most angiosperm genomes, while two small communities (5 and 7) were confined to superasterids and magnoliids, respectively. Of note, taro is the only monocot species with ZAR1 synteny relationship, which was assigned to community 3, along with ZAR1 genes from diverse eudicots, including Caricaceae, Salicaceae, and Rosales (Fig. 6 and Additional Table 7). Next to ZAR1, we also dissected the species constitution of ZRK synteny networks. This analysis pinpointed several relatively conserved communities (1, 5, and 14), covering species from eudicots, monocots, and/or magnoliids. Communities 2, 3, 11, 12, 13, 18, 30, and 36 contained ZRK genes from different genomes belonging to diverse orders of rosids. Moreover, most plant categories have their own specific synteny communities with different sizes. For example, Brassicales and Asterids are represented in multiple communities (5–10), while Amaranthaceae, Rosales, Malpighiales, Malvaceae, and Fagales are only displayed in a few communities (1–3; Additional Fig. 11 and Additional Table 7).

To further zoom into the genomic organization of ZAR1 and ZRK immune gene families and to gain insights into their evolutionary dynamics, we mapped the collinear connections within each of the syntenic communities onto the aforementioned phylogenetic trees (Fig. 7). The results showed strong consistency between synteny characteristics of genomic contexts and phylogenetic relationships of ZAR1 genes. The three ZAR1 groups were found to differ significantly in syntenic conservation. Specifically, group ZAR1-A spanned the most conserved syntenic community (1), while group ZAR1-B was sub-organized by several less interconnected communities (2, 3, and 4). The first splitting branch ZAR1-C showed no syntenic relationships (Fig. 7A). These results suggested that the syntenic context of group ZAR1-B underwent severe rearrangements after splitting from group ZAR1-A. Only a few syntenic connections were identified across the three groups, i.e. between the ZAR1-B group and ZAR1-C from Chimonanthus salicifolius (magnoliids), and between ZAR1-B group and ZAR1-A from taro (monocots), respectively (Fig. 7A). In contrast to the ZAR1-C group, ZAR1-A from magnoliids formed a monophyletic clade interconnected by synteny community 7, indicating the occurrence of ancestral lineage specific transposition events. A similar pattern can be observed in community 5, which consisted of ZAR1-A genes from species belonging to superasterids.

Fig. 7

Phylogenetic classification and syntenic relationships of ZAR1 (A) and ZRK (B) immune genes. The numbers of synteny communities are in accordance with those in Fig. 6 and Additional Fig. 11. The different colored dots on the leaves of phylogenetic trees represent the different orders of the species to which the gene belongs

The relationship between ZRK synteny and phylogenetic classification was found to be much more complex. A considerable proportion of the collinear connections within multiple synteny communities were found to stretch across different ZRK groups (Fig. 7B). These communities with intergroup connections could be caused by background noise as described by Pancaldi et al. [26]. The two largest interrelated communities (1 and 2, each comprising over a hundred connections), are mainly composed of ZRKs from the family Brassicaceae (Fig. 7B, Additional Table 8). The polyploid species B. napus and Camelina sativa are accountable for around 50% of the total collinear relationships among these two clusters, which is larger than the sum of the collinear proportions within diploid species (Additional Table 8). These results indicate that polyploid characteristics contribute largely to the ZRK synteny in Brassicaceae. In addition, syntelogs belonging to several small Brassicaceae-specific communities were found to form monophyletic clades, as well as those from communities confined to Asterids (10 and 26), Rosales (52), Malpighiales (20 and 21), Malvaceae (29), and Fagales (43). These results suggest that ZRKs experienced abundant ancestral transpositions among these lineages.

Discussion

ZAR1 and ZRK homologs across the plant kingdom

Identifying gene homologs can be difficult, especially when analyzing multiple distantly related genomes. Setting a uniform threshold can lead to arbitrary decisions. Here, we combined BLAST similarity searches with phylogenetic analysis to identify ZAR1 and ZRK homologs across the plant kingdom (Additional Fig. 1), focusing on tree topology and phylogenetic relationships rather than relying solely on a BLAST threshold, as used in similar earlier studies [37, 38]. Our method for identifying ZAR1 and ZRKs closely aligns with the approaches used by Adachi et al. [24] and Gong et al. [25]. For those genomes analyzed across these three studies, the numbers of ZAR1-A and/or ZRK genes were largely consistent. However, some minor discrepancies in ZAR1-A copy numbers were observed compared to Adachi et al. (2023), with one additional copy in tobacco, Populus euphratica, and Camelina sativa, and one less in Brassica cretica, cacao, and pomegranate. Similar small differences in copy numbers were found for the ZRKs when compared to Gong et al. (2022), i.e. with one less in taro and Gossypium raimondii, and two less in Solanum tuberosum and Prunus persica. These variations are likely due to different genome versions or updated annotations, as previously noted by others [39].

ZAR1-A–ZRK co-occurrence pattern in early diverging monocot lineages

Taro (Colocasia esculenta) is currently the only monocot species reported to contain both ZAR1 and ZRKs [24, 25]. This study uncovered additional evidence of ZAR1-A–ZRK co-occurrence in other monocots, particularly within early diverging monocot lineages, as demonstrated by the genomes of several Araceae species. Similar as in taro, ZAR1-A and ZRKs were found to co-occur in the closely related aroid species Pistia stratiotes, and this co-occurrence pattern was also observed in Acorus tatarinowii, a species from the first branching monocot lineage Acorales. Amorphophallus konjac, another Araceae member, was found to lack ZAR1-A but to have two ZRK copies. Given that both taro and A. konjac belong to the subfamily Aroideae, we initially hypothesized that ZRKs, rather than ZAR1, might have originated from a shared common ancestor. However, the duckweeds Spirodela polyrhiza and Lemna minor, which belong to the related subfamily Lemnoideae, lack both ZAR1 and ZRKs (Fig. S3, Additional Table 4). This restricted pattern of ZAR1-A and ZRKs co-occurrence to only early diverging monocot lineages (i.e. the alismatid monocots consisting of the orders Acorales and Alismatales) might indicate an ancient interaction or mutual dependence that was preserved throughout their evolution but lost in the lilioids (including Dioscoreales and Asparagales) and commelinids (including Arecales, Zingiberales, and Poales). This observed divergent pattern raises questions not only about the role of complex gene loss events in monocots, but also about whether multiple independent lateral gene transfer events have contributed to the occurrence of ZAR1-A and ZRKs in Araceae. We propose that co-occurrence of ZAR1-A and ZRKs can also be anticipated in species belonging to other early monocot clades.

De-coevolution of the ZAR1-A–ZRK co-occurrence pattern in Ericales

Based on their observation of consistent co-presence or absence patterns of ZAR1-A and ZRK genes across the plant genomes examined, Gong and colleagues concluded that the ZAR1-ZRK resistosome originates before the divergence of eudicots and monocots. This co-occurrence pattern of ZAR1-A and ZRKs was confirmed in most of the genomes studied in our work, but with clear exceptions on co-presence/co-absence in several genomes belonging to Ericales that were found to contain ZRKs but lack ZAR1-A (Additional Table 4). Ruling out that the absence of ZAR1-A in the five genomes of Ericales (Vitellaria paradoxa, Camellia sinensis, Actinidia chinensis, Vaccinium darrowii, and Rhododendron simsii) is due to low-quality of genome assembly (BUSCO values ranging from 85.3 to 91.0%), this would suggest de-coevolution of the ZAR1-A–ZRK co-occurrence pattern in angiosperms. An opposite pattern in which species contain ZAR1-A but lack ZRKs was, however, not found in this study (Additional Table 4). These observations raise several new questions. Which roles do ZRKs play in species that lack ZAR1-A, and are these ZRKs recruited by other NLRs to form immune complexes, or do they play functional roles in distinct immune mechanisms? We hypothesize that in Ericales species lacking ZAR1-A, ZRKs may act as effector sensors in resistosome interactions with NLRs other than ZAR1 or function independently of NLRs as decoys for kinase virulence targets to inhibit infection [40]. This could explain why these species retain ZRKs instead of losing them entirely. The role of these ZRKs, when not associated with ZAR1, can potentially be elucidated through comparative functional and biochemical analysis [23,24,25].

Potential factors influencing ZAR1-A–ZRK co-occurrence patterns

The variations on the ZAR1-A–ZRK co-occurrence pattern observed across several plant lineages could be attributed to multiple factors. Next to gene loss that could affect copy number, genes may have functionally diverged, resulting in altered interactions or regulatory networks across different lineages. Such functional divergence might result from adaptation to specific ecological niches or shifts in environmental pressures. Emerging perspectives are shedding new light on the fate of genes retained as duplicates after polyploidization. For example, a recent study revealed functional divergence in CIPK6 homologs in upland cotton, where GhCIPK6D1 promotes drought sensitivity, while GhCIPK6D3 enhances drought tolerance [41]. Gene mobility might also have had an impact, with horizontal gene transfer or mobile genetic elements potentially affecting the presence and copy number of ZAR1 and ZRKs in certain lineages, leading to discrepancies in their co-occurrence pattern.

Conclusion

In this study, we explored the phylogenomic relationships of ZAR1 and ZRK immune gene families using an extensive set of plant genomes—covering the large diversity in plant lineages. Our results confirmed an early origin and a co-evolution pattern of ZAR1-A and ZRKs. Compared to previous studies, our study revealed novel ZAR1-A–ZRK co-occurrence in monocots, particularly in the early diverging lineages Araceae and Acorales. Next, we found several robust exceptions to the co-evolutionary ZAR1-A–ZRK pattern in genomes of Ericales, implying de-coevolution driven by divergence from the shared co-evolution mode. Our study also explored the syntenic contexts of these two immune gene families, revealing both conservation and lineage-specific patterns. These findings expand our knowledge on the evolutionary history of the ZAR1 and ZRK immune system.

Methods

Genome searches and gene identification

A total of 231 plant genomes was used in this study, of which most were previously used by Pancaldi et al. [26] (Additional Table 1). Taxonomic classification of species was obtained using APG website and NCBI database [1]. The species taxonomy tree was constructed using ETE 3.1.1 followed by visualization in iTOL v5 [42]. Amino acid sequences of ZAR1 genes from Arabidopsis thaliana, Nicotiana benthamiana, Solanum lycopersicum, and Manihot esculenta [23], and 13 ZRKs from Arabidopsis thaliana were used as queries for reciprocal BLAST 2.14.0 searches against the protein sequence databases with an e-value threshold of 1e-2 [14]. Obtained amino acid sequences were subsequently compiled and aligned using MAFFT v7 webtool [43]. The resulting multiple sequence alignment (MSA) was filtered by trimAl 1.2rev59 to remove gapped columns (parameters: gt 0.8, st 0.001, cons 60), followed by maximum likelihood tree construction in FastTree 2.1.11 [44]. Only hits clustering with query sequences and having a clear evolutionary relationship were retained and designated as true homologs. Details on the bioinformatic pipeline can be found in Additional Fig. 1.

Phylogenetic analysis and gene classification

Amino acid sequences of 254 ZAR1 and 929 ZRK genes were used for phylogenetic analysis and classification. The phylogenetic tree of the ZAR1 family was rooted using three homologs (XP_006829323.1, XP_006836665.2, and XP_020520999.1) from Amborella trichopoda, a species representing the first-branching lineage of angiosperms. The phylogenetic tree constructed of ZRK genes was rooted using ATA10.962.t1 from Acorus tatarinowii, a species belonging to the earliest branching lineage of monocots. MSA was generated by MAFFT v7 using default settings and filtered using trimAl 1.2rev59 to remove gapped columns as described above (Additional File 1–2) [43, 45]. Subsequently, maximum-likelihood tree was generated by IQ-TREE v2.2.2.9 with the settings of MFP model and 1000 bootstrap replicates [46], and finally visualized in iTOL v5 [42].

Characterization and chromosomal distribution of gene duplications

Chromosomal locations of ZAR1 and ZRK genes were extracted from the corresponding genome annotation gff3 files by an inhouse Perl script and visualized in Mapchart software [47]. The duplicate_gene_classifier function integrated in MCScanX was used to assess gene duplication patterns, including singleton, dispersed, proximal, tandem, and WGD or segmental duplication [48]. Gene pairs derived from WGD or segmental duplication in B. napus were visualized using Circos software [49].

Selection pressure analysis

Ratios of nonsynonymous substitutions (Ka) to synonymous substitutions (Ks) of homologous pairs were calculated using MUSCLE and KaKs_Calculator implemented in ParaAT2.0 software [50,51,52]. Ka/Ks values < 1 indicate negative or purifying selection, while Ka/Ks values > 1 represent positive selection.

Synteny network analysis

All-against-all comparisons of protein sequences were performed using Diamond v2.0.11.149 software [53]. The top five blast hits were retained as inputs for the detection of syntenic blocks using MCScanX (minimum match size 3; maximum gaps 25) [48]. Outputs were organized in a synteny network, with nodes representing genes and edges representing collinear relationships between genes. The edges with both ZAR1/ZRK gene nodes were extracted from the whole synteny database and represented as the ZAR1/ZRK synteny network, which was finally visualized in Gephi v0.9.2 [54]. Infomap implemented in the R package igraph was used to identify syntenic communities [55], and the retrieved communities consisting at least three nodes were graphically represented in Cytoscape v3.8.2 [56]. Additionally, species and gene composition within the communities were dissected using an inhouse R script. Based on gene copy numbers per species, Jaccard method was used to calculate the dissimilarity indexes of all communities, and ward.D was employed for hierarchical clustering [57, 58]. Syntenic connections between nodes of identified communities were graphically displayed in phylogenetic trees using iTOL v5 [42].