Diversity and Evolution of NLR Genes in Citrus Species

1. Introduction

When pathogens infect plants, the plants deploy two innate immune systems: pathogen-associated molecular pattern-triggered immunity (PTI) and effector-triggered immunity (ETI), to defend against the invasion of pathogens [1,2]. PTI serves as the primary line of defense, promptly recognizing and responding to a wide range of pathogen-associated molecular patterns (PAMPs) or microbe-associated molecular patterns (MAMPs). This triggers fundamental immune mechanisms such as ethylene production, oxidative bursts, callose deposition, and the expression of defense-related proteins [3,4], thereby impeding the invasion and spread of pathogens. However, pathogens evolve various strategies to counteract or evade PTI defenses. One common strategy is the release of effector molecules [5], which can interfere with or inhibit key components of PTI, aiding pathogens in colonizing and proliferating within the plant host.

Fortunately, plants possess a second, more specialized immune system: ETI [5]. When pathogens release effectors to suppress PTI, the plant’s resistance genes (R genes) can directly or indirectly recognize these effectors, activating ETI immune responses. ETI typically leads to a hypersensitive response (HR) at the infection site, where cells rapidly undergo programmed cell death, limiting the further spread of pathogens [6]. These two immune systems are interconnected, working together to maintain the health and survival of plants [2,7,8].

R genes are critical components of the immune system’s ETI, playing a key role [6]. Most R genes encode disease resistance proteins with nucleotide-binding site and leucine-rich repeat (NBS-LRR) domains—these genes are thus often termed NLR genes [9]. The R genes discussed in this paper are NLR genes, excluding other types. Beyond the canonical NBS and LRR domains, the N-termini of R proteins frequently feature additional domains or motifs, such as the Toll/Interleukin-1 receptor (TIR) domains, the powdery mildew resistance protein 8 (RPW8) domains, and coiled-coil (CC) motifs [1,10]. Proteins harboring the TIR and RPW8 domains are classified as TNL and RNL, respectively [10,11]. Similarly, proteins with CC motifs are designated as CNL [9]. A significant proportion of R proteins, possessing only the NBS and LRR domains, are collectively categorized as NL (NBS-LRR).

Comparative genomics research has shown that NLR genes are extensively distributed among land plants and charophytes, indicating a deep evolutionary origin [1,10,12]. Additionally, NLR genes have been identified in a limited number of green algae species, further supporting their ancient lineage [13]. Phylogenetic analyses suggest that the emergence of R genes preceded the evolutionary split between charophytes and land plants [1,10]. Following their origin, NLR genes have undergone adaptive evolution, significantly contributing to the transition of plants from aquatic to land habitats [10].

Citrus is one of the world’s most significant fruits, providing essential nutrients such as vitamin C and citric acid. Cultivation of citrus spans over 140 countries and regions globally, thriving in warm and moist climates below the northern latitudes [14]. In 2020, the global citrus cultivation area reached 521,000 hectares, yielding a production of 96.81 million metric tons [15]. However, citrus faces various threats from diseases, such as citrus HLB, a devastating ailment that severely impacts the yield and quality of citrus. Addressing this challenge necessitates a deep understanding of citrus’s mechanisms of disease resistance. The study of the diversity of disease-resistant genes is foundational for elucidating these defense mechanisms. Currently, there is no systematic research reported on citrus NLR genes.

This study integrates comparative genomics and phylogenetics to systematically explore the NLR genes in the citrus genome, unraveling the mechanisms underlying the diversity, origin, and evolutionary patterns of citrus NLR genes. We observed substantial diversity in NLR genes among different citrus species. Further in-depth investigation revealed that gene duplication and recombination are the primary mechanisms shaping the diversity of NLR genes in citrus. This study holds significant implications for a profound understanding of citrus NLR genes and for breeding disease-resistant citrus varieties.

2. Materials and Methods

2.1. Data Used in This Study

At the outset of this study (September 2020), genomic data for 10 citrus species were accessible, namely C. medica, C. ichangensis, C. hindsii, Atlantia buxifolia, C. trifoliata, C. grandis, C. mangshanensis, C. sinensis Osbeck cv. newhall, C. sinensis Osbeck cv. Valencia, and C. clementina. The genome of C. sinensis Osbeck cv. newhall was assembled and annotated by our research team [29], available for download from Zenodo database (https://zenodo.org/record/8176803). Genomic data for the remaining 9 citrus species can be obtained from the CPBD (http://citrus.hzau.edu.cn/download.php) [22].

2.2. Identification of NLR Genes

The nine seed sequences of the NB-ARC domain (PF00931) were obtained from the European Bioinformatics Institute (EBI) database and stored in Stockholm format. A Hidden Markov Model (HMM) of NB-ARC was constructed using HMMER package 3.1b2 [30]. Then, this model was used to search for proteins containing NB-ARC domain among the protein sequences of 10 citrus genomes (Table S1). The Blastp algorithm was also employed to search for citrus proteins containing NB-ARC domain, using the NLR proteins of Arabidopsis thaliana [10] as queries, with a cutoff E-value of 10^-3. The results obtained from both methods were merged and redundant sequences were removed to obtain citrus proteins containing NB-ARC domain. Multiple sequence alignment of citrus proteins containing the NB-ARC domain and previously published NB-ARC or NAT domain proteins from both eukaryotes and prokaryotes [31] was performed using MAFFT software [32] with G-INS-i algorithm, followed by manual edited using MEGA11 [33]. ModelFinder [34] was utilized to search for the best model, and IQ-TREE software [35] was employed with JTT+R10 model for phylogenetic analysis. Proteins with MalT and SWACOS domains were used as outgroups [31]. The protein domain architecture was annotated using CD-Search [36]. NLR genes were searched in the nucleotide sequences of 10 citrus genomes using NLR-Annotator software [37] with default parameters, for testing analysis of the identification results.

2.3. Phylogenetic Analysis

To investigate the origin of NLR genes in citrus, 3168 NLR gene sequences from 16 species including citrus, Sapindaceae, and Anacardiaceae, were aligned using MAFFT [32] with G-INS-i algorithm, manual edited using MEGA11 [33]. The ModelFinder [34] was used to research the best model, and then Q.mammal+R10 was employed for phylogenetic tree construction using IQ-TREE [35].

The citrus NLR gene sequences were aligned using MAFFT [32] for multiple sequence alignment, followed by manual edited using MEGA11 [33]. The ModelFinde [34] was employed to identify the best model, and subsequently, JTT+F+R10 was used for phylogenetic tree construction using IQ-TREE [35]. The phylogenetic analysis of citrus species based on whole genomes was conducted using the STAG algorithm implemented in OrthoFinder [19].

2.4. Positive Selection Analysis

Due to the significant divergence of NLR genes, alignment and positive selection analysis were not feasible for all sequences. Consequently, we divided them into 65 groups based on phylogenetic analysis (Figure 4). Nucleotide sequences were aligned using MAFFT [32] and manually refined with MEGA11 [33]. A likelihood ratio test for positive selection was conducted by comparing the M7 (beta) and M8 (beta&v) models using the PAML4.9 software package [38,39]. For gene groups demonstrating positive selection via the likelihood ratio test, positively selected sites were identified through the Bayes empirical Bayes (BEB) procedure.

2.5. Recombination Analysis

Recombination was detected using the RDP5 package [40]. The aligned nucleotide sequences were imported into RDP5 and analyzed with the RDP, GENECONV, Chimaera, MaxChi, BootScan, SiScan, and 3Seq methods. Potential recombination events were deemed significant only if supported by at least four methods with a p-value < 10⁻⁶.

3. Results

3.1. Identification and Diversity of Citrus NLR Genes

This study utilized comparative genomics and phylogenetics to identify NLR genes in the proteomes of 10 citrus species. Initially, protein sequences were screened for the NB-ARC domain using the HMM model, resulting in the identification of 1,875 sequences containing the NB-ARC domain (Table S1). Subsequently, 1,845 sequences with the NB-ARC domain were identified via Blastp homology searches (Table S1). After redundancy reduction, sequences containing the NB-ARC or NAT domains were used to construct a phylogenetic tree, along with reported homologous sequences from representative plants and other species (including animals, bacteria, fungi, and archaea). We identified citrus proteins with the NB-ARC domain that formed monophyletic groups with plant R genes, designating these as citrus NLR protein sequences (Figure 1A). Some sequences with the NB-ARC domain formed monophyletic groups with non-plant homologs and were therefore excluded. Ultimately, we identified 1,585 NLR genes (Table S1).

Utilizing the HMM method and the Pfam database, we annotated the domain structures of the 1,585 citrus NLR genes, identifying 183 as TNL and 35 as RNL (Figure 1B). Notably, we also annotated the Rx_N (PF18052) domain, predicted to be a CC motif but suggested by research to be a four-helical bundle fold [16]. The Rx_N domain is equivalent to the CC domain annotated in previous studies [16]. Thus, in this context, RxN-NBS-LRR corresponds to CNL, as identified in prior research [9,17]. We collectively termed the NLR proteins containing the Rx_N domain as CNL, totaling 272 proteins (Figure 1B). Additionally, a substantial number of NL proteins containing only NBS or NBS-LRR domains were annotated, totaling 1,095 (Figure 1B).

Among the 10 citrus species, aside from the 108 TNL sequences in C. clementina, the number of TNL sequences in the other nine species was relatively low. Only a few were found in C. medica, C. ichangensis, C. hindsii, and C. trifoliata, while they were absent in C. mangshanensis and C. grandis (Figure 1B). RNL was the least abundant among the three types, evenly distributed across all species, ranging from 3 to 4 sequences each (Figure 1B). There were six RNL sequences in C. sinensis Osbeck cv. Newhall, possibly due to its diploid nature, resulting in a doubling of the quantity.

3.2. The Origin of Citrus NLR Genes

The study utilized 1,585 identified NLR proteins and 1,583 published NLR proteins from the Anacardiaceae and Sapindaceae families (Table S2) to construct phylogenetic relationships (Figure 2). The results indicate that citrus NLR proteins fall within the diversity of Sapindaceae, suggesting that citrus NLR genes originated comparatively later than those in Sapindaceae. Within the Sapindales order, Rutaceae is the oldest, followed by Sapindaceae, with Aurantioideae (Citrus genus) being the most recent [18]. In most monophyletic groups, the phylogenetic relationships of NLR proteins align with the species relationships of Rutaceae, Sapindaceae, and the Citrus genus (Figure 2). The coevolutionary relationships between NLR proteins and species suggest that NLR genes originated concurrently with the species’ origins.

3.3. The evolution of citrus NLR genes

We employed maximum likelihood methods to construct phylogenetic relationships among citrus NLR proteins. Integrating domain architecture of NLRs, we identified five main clusters in the phylogenetic tree of citrus NLR proteins, comprising one cluster of TNLs, one cluster of CNLs, one cluster of RNLs, and two clusters of NLs (Figure 3). NLs exhibited the highest diversity, while TNLs, RNLs, and CNLs fell within the diversity of NLs, indicating their origin from NLs with the acquisition of TIR, RPW8 domains, and CC motifs respectively. Additionally, NLs were found within both TNLs and CNLs clusters (Figure 3), indicating TNLs and CNLs lost their TIR and Rx_N domains, respectively, transforming into NLs. Overall, citrus NLR genes underwent frequent events of domain acquisition and loss during evolution.

To test whether citrus NLR genes underwent adaptive evolution after their origin, positive selection analysis was conducted on 65 groups selected base on phylogenetic tree. amino acid residues under positive selection was detected in the majority of groups (52/65)(Figure 4). Among the 52 gene groups where positive selection was detected, the number of positively selected sites varied(Table S3). Our findings indicate that the vast majority of citrus NLR clusters underwent adaptive divergence.

3.4. Mechanism of NLR Genes Diversity Formation in Citrus

The phylogenetic tree of citrus was constructed using OrthoFinder software [19] based on full genome of 10 citrus genomes (Figure 1A). Atlantia buxifolia is placed at the root of the tree, with the overall phylogenetic order largely consistent with previous reports [20]. Comparing the number of NLR genes across different citrus species on the phylogenetic tree reveals considerable diversity in NLR gene counts even within the same subfamily. Among the 10 citrus species, substantial differences exist in the copy numbers of NLR genes. For instance, C. clementina has 405 copies, whereas C. ichangensis has only 56, making the former about eight times more abundant than the latter. In the case of the sweet orange subspecies C. sinensis Osbeck cv. Newhall and C. sinensis Osbeck cv. Valencia, 178 and 107 copies were identified, respectively. Notably, while C. sinensis Osbeck cv. newhall is a diploid genome, C. sinensis Osbeck cv. Valencia is a haploid genome. Although the genomic content of C. sinensis Osbeck cv. Newhall is double that of C. sinensis Osbeck cv. Valencia, the NLR gene copies are not proportionally higher. What evolutionary mechanisms have contributed to the genetic diversity of citrus NLR genes?

Based on the phylogenetic tree of citrus NLR proteins, we investigated horizontal gene transfer and gene duplication events within citrus NLR genes. Gene duplication events were identified when NLR genes from the same species formed a monophyletic group on the phylogenetic tree. Horizontal gene transfer events were identified when the genetic diversity of one species fell within that of another species.

Following statistical analysis, we generated a diagram depicting the levels of horizontal gene transfer (HGT) between various citrus species (Figure 5A), revealing relatively frequent gene exchange among C. sinensis Osbeck cv. Newhall, C. sinensis Osbeck cv. Valencia, C. mangshanensis, and C. clementina. Notably, we observed gene exchange between C. sinensis Osbeck cv. newhall and C. clementina with Atlantia buxifolia. Additionally, there was minimal gene exchange between C. sinensis Osbeck cv. Newhall and C. medica, while the remaining four citrus species (C. ichangensis, C. hindsii, C. trifoliata, and C. grandis) and other species showed no NLR gene exchange. Correlation analysis between the number of NLR genes and the frequency of horizontal gene transfer events revealed a non-significant positive relationship (Figure 5B), suggesting that horizontal gene transfer contributes insignificantly to citrus NLR gene diversity.

Recombination has long been considered a significant driver for NLR gene diversification [21]. We used seven algorithms within the RDP software (Martin et al., 2015) to search for recombination signals within citrus R genes. A single gene recombination event was considered when four or more algorithms detected positive signals. In this study, 173 recombination events were identified, with C. clementina exhibiting the highest number at 81 occurrences (Table S4). Correlation analysis between the number of NLR genes and the number of recombination events across citrus species showed a significant positive relationship (Figure 5C), indicating that gene recombination is one of the mechanisms contributing to citrus NLR gene diversity.

A total of 293 gene duplication events were recorded in this study, with occurrences ranging from 6 to 116 across different citrus species. C. clementina once again exhibited the highest number of gene duplication events at 116 occurrences (Table S4). Correlation analysis between the number of NLR genes and the number of gene duplication events revealed a positive relationship, indicating that gene duplication also contributes to citrus NLR gene diversity. In summary, gene recombination and gene duplication are the main mechanisms shaping NLR gene diversity in citrus.

4. Discussion

In this study, we employed comparative genomics and phylogenetic methods to identify NLR genes across 10 citrus species. To validate the effectiveness of our mining approach, we used the NLR-Annotator software (Zhang, 2020) to search nucleotide genomes, resulting in the retrieval of 1847 citrus NLR gene sequences. The NLR counts extracted by NLR-Annotator across various citrus species closely matched the results obtained from Blastp and HMM of this study. Comparison with previously reported citrus NLR gene numbers in the ANNA database (Liu et al., 2021) revealed close correspondence between the quantities in A. buxifolia, C. grandis, C. medica, and C. ichangensis and the results obtained in this study through Blastp and HMM. Thus, the identification methods utilizing Blastp and HMM employed in this study are validated.

The ANNA database reports 529 NLR genes of C. sinensis Osbeck cv. Valencia [23], whereas our study’s Blastp and HMM analyses yielded 179 and 178, respectively, and the NLR-Annotator software result was 181. The substantial variance can be attributed to the utilization of the C. sinensis Osbeck cv. Valencia v2.0 genome [24] in prior studies, while we employed the v3.0 genome version [25]. Factors such as genome assembly quality and annotation precision between the two genome versions may influence the NLR gene number. We conducted NLR gene mining on the nucleotide sequences of both sweet orange genome versions using the tblastn algorithm. We found comparable outcomes with no significant differences between the two genome versions. This underscores annotation precision as the primary factor influencing the quantity of NLR genes.

We integrated the results from Blastp and HMM, eliminating redundancies. Subsequently, we employed phylogenetic analysis to further screen the sequences, remove false positives, and discard NLR sequences generated by alternative splicing, thereby obtaining the NLRs for each citrus species in this study. We employed a more stringent identification method to ensure the accuracy of the results.

The distribution of TNLs varies significantly among citrus species. TNLs are sparsely present or absent in C. medica, C. ichangensis, C. hindsii, C. trifoliata, C. grandis, and C. mangshanensis, while in C. clementina, their number can be as high as 108. To exclude the impact of gene annotation on the results, we used the tblastn algorithm to mine the genomes of C. medica, C. ichangensis, C. hindsii, C. trifoliata, C. grandis, and C. mangshanensis for TNLs. The identification results were consistent with our study, thus eliminating the influence of genome annotation on the identification results. The evolutionary fate of NLR genes is mainly determined by natural selection and genetic drift [21]. In this study, positively selection amino acid residues were detected in all four clusters of TNLs selected, suggesting that the loss of TNLs in species such as C. mangshanensis may be due to genetic drift.

In-depth analysis revealed that the expansion of TNL gene copies in C. clementina is due to frequent gene recombination and duplication events. On one hand, previous studies have speculated multiple origins of TNLs [26], while loss has been reported in Aquilegia coerulea and monocots [27]. On the other hand, research has indicated that the evolutionary rate of TNLs is higher than that of non-TNLs [26], suggesting that, compared to CNLs and RNLs, TNLs may be more prone to acquiring or losing the TIR domain, thus undergoing a more complex evolutionary process.

Differences in NLR gene copy numbers among citrus species are not uncommon, as seen in the Rosaceae family, where the wild strawberry (Fragaria vesca) has 144 copies, while the apple (Malus pumila) has 748 copies [26]. This study further analyzes gene recombination and duplication as the primary mechanisms shaping NLR gene diversity in citrus. While horizontal gene transfer did not exhibit significant impact across all citrus species, it was the major driver of NLR gene copy number increase in A. buxifolia. The number of gene duplications and recombinations in C. medica and A. buxifolia is the same, with 10 and 6 instances, respectively. Horizontal gene transfer occurred 9 times in A. buxifolia and only once in C. medica. The NLR gene copy number in A. buxifolia is 1.67 times that in C. medica. Therefore, horizontal gene transfer is the primary mechanism driving the increase in NLR gene copy number in A. buxifolia.

Approximately 30 million years ago, the Rutaceae subfamily originated from the ancient Indian plate, which was then connected to Africa and Australia. As the plates collided and drifted, the Rutaceae subfamily gradually diversified, with the citrus genus originating over 90 million years ago in southern China [20]. During the colonization of terrestrial environments, plants encountered new pathogen species and infection levels [28]. Positively selected NLR genes play a crucial role in plant-pathogen interactions. This evolutionary process involved frequent gains and losses of TIR, RPW8, and LRR domains, enhancing R-protein diversity and enabling adaptation to various environmental challenges [10]. The discovery of numerous adaptively evolved NLR genes in citrus suggests their significant role in global citrus colonization.

5. Conclusions

This study provides a comprehensive understanding of the diversity and evolution of NLR genes in citrus species, which were previously underexplored compared to other plants. By identifying 1,875 NLR genes across 10 citrus species and analyzing their variations, origins, and adaptive evolution. The findings are particularly significant for their potential application in breeding disease-resistant citrus varieties, making the study a foundational resource for future citrus genetic research and crop improvement efforts.