In this study, we selected three Australian limes and two common citrus cultivars ( Table 1 ). In the Australian limes, C. glauca was considered as HLB-resistant and in the HLB-resistance category C2, according to Ramadugu’s (Ramadugu et al., 2016b) evaluation system, whereas C. australasica and C. inodora were tolerant, in category C3. The two citrus cultivars C. clementina and C. sinensis were both susceptible and rated as C7. The chromosomal-scale genome assemblies of the three Australian limes are available at the NCBI Sequence Read Archive (SRA) (Singh et al., 2024), with their GenBank assembly numbers being GCA_029618585.1 (C. australasica), GCA_029633175.1 (C. glauca) and GCA_029721495.1 (C. inodora). The genome sequences of C. clementina and C. sinensis were downloaded from the Citrus Pan-genome to Breeding Database (http://citrus.hzau.edu.cn/orange/) (Wang et al., 2021; Wu et al., 2014), and their GenBank assembly numbers at NCBI are GCA_000493195.1 and GCA_019144185.1, respectively. These assemblies were of similar genomic sizes (298.9-376.5 Mb) and assembly quality with the scaffold N50 ranging from 28.9-37 Mb and scaffold L50 being 4 to 5. To ensure consistency of chromosome numbering between genomes, genomic alignment between C. sinensis and C. clementina was inspected using a web-based genome compare tool D-GENIES (Cabanettes and Klopp, 2018) (https://dgenies.toulouse.inra.fr). Chromosome numbers on C. sinensis were reassigned (for the present study) according to C. clementina, which was used as reference to order the scaffolds of the three Australian lime genome assemblies (Singh et al., 2024).

The primary haplotype assembly of the three Australian limes for each genome sequence was processed to remove any potential soft masking using in-house awk script (lowercase ACTGs that signify an annotated repeat or transposable elements). We downloaded the reference RefPlantNLR database (Kourelis et al., 2021) and the meme.xml file from the NLR-Annotator v.2 package (Bailey et al., 2015; Steuernagel et al., 2020). We executed the FindPlantNLR snakemake pipeline (Chen et al., 2023) to identify and annotate the NLR genes based on 13 Pfam accessions ( Supplementary Table S1 ). The FindPlantNLR pipeline relies on tblastn, nhmmer, and NLR-Annotator to identify NLR loci (Camacho et al., 2009; Wheeler and Eddy, 2013; Steuernagel et al., 2020). The pipeline extended the potential region and 20 kb flanking sequence with bedtools (Quinlan and Hall, 2010). The extracted sequence was then input into the BRAKER annotation pipeline using the RefPlantNLR database as a model (Hoff et al., 2016). The resulting genes were scanned for motifs using Interproscan (Jones et al., 2014) followed by classification and script annotation within the FindPlantNLR package.

To determine which FindPlantNLR generated R-genes were retrieved from unannotated regions in the original genome annotation, we compared the NBARC annotation (.gff file) and the corresponding genome annotation using the gffcompare (v.0.9.12) software package within the GFF utilities (Pertea and Pertea, 2020). The results of gffcompare reported the matched and unmatched transcripts between the two annotation files. The R-genes that found no match in the genome annotation were considered as masked R-genes, i.e. they were not categorized as genes in the original genome annotation. The matched genes were categorized based on the types of matching relationship to reference transcripts and indexed by classification codes (detailed in Supplementary Table S2 ).

To compare the complement of R-genes identified in this study with previously reported R-genes in C. sinensis and C. clementina, we first made BLAST databases from the protein FASTA files of this study using makeblastdb. We compared the protein sequences of R-genes using blastp with HSP (high scoring pair) e-value set at 1e^-6.

Phylogenetic analysis of the five citrus species was performed on the protein sequences of the NBARC genes using the multiple sequence alignment program MAFFT (v7.526) (https://mafft.cbrc.jp/alignment) (Katoh and Standley, 2013) with the default setting. The phylogenetic relationship was visualized in a circular plot using the visualization program Chiplot (https://www.chiplot.online) (Xie et al., 2023).

For each citrus species genome, the BUSCO (v5.8.1) (Benchmarking Universal Single-Copy Orthologs) genes were retrieved and concatenated to generate an alignment supermatrix using a Python pipeline (https://github.com/jamiemcg/BUSCO_phylogenomics), which was used to construct species phylogeny using IQ-TREE (v2.3.6) (Nguyen et al., 2015) with 1000 bootstrap replicates. The tree structure was visualized using MEGA11 (Tamura et al., 2021).

Syntenic genes between each pair of citrus species were identified using the One StepMCScanX-SuperFast module in TBtools (v2.086) (Chen et al., 2020). Among these syntenic genes, 39 were identified to be common across all species and were used to construct phylogenetic trees following alignment and concatenation. The distribution and organization of these consensus syntenic genes were visualized using Multiple Synteny Plot in TBtools (Chen et al., 2020). Phylogeny based on these consensus syntenic genes was constructed using IQ-TREE (Nguyen et al., 2015) after alignment and concatenation. The syntenic genes within the genome of each species were analyzed in a similar way as described above and were visualized using the Advanced Circos module in TBtools (Chen et al., 2020). The clusters were defined by the presence of at least three genes that were located less than 200 kb apart (Holub, 2001).

To estimate the selection pressure acting on R-genes, we calculated the rates of nonsynonymous substitution, synonymous substitution, and their rate ratio (Ka, Ks, and Ka/Ks) on each pair of syntenic R-genes between species using DnaSP6 (Rozas et al., 2017). Loci with calculated Ka/Ks values were plotted on chromosomes using the R package chromoMap (v0.4.1) (Anand and Rodriguez Lopez, 2022).

The protein sequences were employed to infer structural and functional differences in R-genes between genomes using Orthovenn3 (Sun et al., 2023) with the OrthoMCL algorithm (E-value = 1e-2, inflation value = 1.5). Pairwise sequence similarities between species were calculated with a threshold of the e-value ≤ 1e⁻⁵, and the inflation value was set at 1.5 for orthologous cluster generation using the Markov clustering algorithm. The expansion or contraction in gene family sizes was analyzed using the CAFE5 (Mendes et al., 2021), which implements a birth-death model to infer phylogenetic history and evolutionary time. The biological processes and molecular functions associated with identified gene clusters were retrieved and identified with GO terms annotation. Based on the identified unique clusters, the unique genes were collected and used for GO enrichment analyses using the built-in feature in OrthoVenn3.

The sequences and coordinates of R-gene domains, including CC, TIR, and LRR, were retrieved using Hidden Markov model (HMM). Briefly, the HMM profiles using HMMER (v3.4) (http://hmmer.org/) were mapped on protein sequences (Singh et al., 2024) from C. australasica, C. glauca, C. inodora, C. clementina and C. sinensis. Accessions for CC, TIR, and LRR were retrieved from PFAM database version 36.0 (http://pfam.xfam.org/). Six LRR domains were successfully retrieved and designated as: LRR1 (PF00560.37), LRR3 (PF07725.16), LRR4 (PF12799.11), LRR5 (PF13306.10), LRR6 (PF13516.10), and LRR8 (PF13855.10). To ease the visualization of the domain organizations on R-genes, we curated 20 genes from each species that best represent the complete gene set of each species. In this process, we first retrieved and combined the sequences that contain all the R-gene domains using an in-house program (https://github.com/saikizu/DoBioPython), followed by alignment using Clustal Omega (v1.2.4) (Sievers et al., 2011). The alignments were used to construct phylogenetic trees using IQ-TREE (v2.3.4) (Minh et al., 2020) and the tree files (in Newick format) were then used in the python-based program Treemer (Menardo et al., 2018) to prune the leaf numbers down to 20 by eliminating those that contribute the least to the tree diversity. The domain composition and organization were visualized using Simple Biosequence Viewer in TBtools, v1.108 (Chen et al., 2020).

EMBOSS-CONS (v 6.6.0.0) (Madeira et al., 2024) was used to generate consensus sequences from the domain sequences in each species to compare the domains between species. In this process, a residue is considered to be consensus if the number of positive matches at the position is greater than half (≥3) of all the sequences in the alignment.

Motif detection was performed using the MEME SUITE (v5.5.7) (Bailey et al., 2015) using protein sequences of the NLR genes. Each sequence identified the top five motifs with a motif width between 6 and 50 amino acids. Identified motifs were concatenated, aligned, and subjected to phylogeny construction for comparison.

The D-GENIES (Cabanettes and Klopp, 2018) alignment was used to reassign chromosome numbers to genomes that were labelled differently to facilitate comparison between species. In the de novo assembled genomes of the three Australian limes, C. clementina was used as the reference for numbering the chromosomes. Mapping between C. sinensis and C. clementina identified 69.61% of matches with greater than 50% similarity between the two genomes ( Supplementary Figure S1A ). The dot plots ( Supplementary Figures S1B, C ) indicated five chromosomes (Chr1, 3, 4, 5, and 7) in C. sinensis were inconsistent with C. clementina and were relabeled (for the present study) as follows: Chr1→Chr7, Chr3→Chr5, Chr4→Chr1, Chr5→Chr3, Chr7→Chr4.

Using the R-gene specific annotation pipeline FindPlantNLR (Chen et al., 2023), we identified a wide range of NBARC genes in the five citrus species ( Table 1 ). These NBARC genes accounted for about 1.4-2.5% of the total genes in each genome. The two cultivars harbored 11.8-88.4% more NBARC genes than each of the Australian limes. The lowest number of NBARC genes was found in C. inodora. Comparing the R-gene annotation with the original genome annotation, we found that FindPlantNLR uncovered many NBARC genes that were not predicted in the original gene annotation of each genome. Over 70% NBARC genes were not previously annotated in C. australasica and C. glauca and 30-40% had not been annotated in the other three species. For the NBARC genes that showed overlaps in the genome annotation, the matching relationship were categorized into 11 classes, in which three (coded as j, k, and o) represented more than 50% of all the NBARC genes ( Supplementary Table S2 ). To evaluate the completeness of R-genes predicted in C. clementina and C. sinensis in this study, we compared the protein sequences of these newly annotated R-genes with the R-genes annotated in previous studies (Wang et al., 2015; Yin et al., 2023). The BLAST results indicated that the R-genes identified in this study contained the complements of previously reported R-genes in both species ( Supplementary Tables S3 - S5 ), with the average identity more than 94.8%.

Among the NBARC genes, most (78.9-86.2%) were found to harbor an LRR as well, thereby classified as NLR genes. Similar to NBARC genes, the two cultivars had 15.8-104% more NLR genes than each of the three Australian species, and, C. inodora had the fewest. About 50% of NLR genes in each species fell into the CNL or TNL category, both of which varied highly between species. While the numbers were similar in C. glauca, the CNL genes were nearly two times higher than that of TNL genes in C. clementina and 47.5-63.8% higher in other species. With a few NLR genes identified as RNL, the rest of the NLR genes were classified as NL, which accounted for about 50% of the total NLR genes.

Phylogenetic analysis using protein sequences of the NBARC genes of the five citrus species generated three major clades ( Figure 1A ), which corresponded to three major R-gene families, namely CNL (red), TNL (blue) and NL (green). Other NBARC type of genes (orange) formed small clusters and dispersed among the phylogenic branches. In contrast to the distinct clades between R-gene classes, the relatedness of the species were not adequately depicted, and there was no distinction between the three Australian species, C. australasica (purple), C. inodora (orange) and C. glauca (red), and the two cultivated species of C. clementina (pink) and C. sinensis (green).

Using the protein sequences, we first explored the relationship between the five citrus species on a whole genome scale. We identified and aligned all single-copy BUSCO genes in all the genomes and constructed a robust maximum likelihood phylogenetic tree ( Figure 1B ). The phylogeny grouped the five species into two major clades, with C. australasica and C. inodora forming one clade and the other three species forming the other clade, in which the two cultivars form a sub-clade. It’s noticeable the support values on most branches were low, signifying a lack of strong distinction between the clades represented by the BUSCO genes.

To identify the conservation of homologous genes and R gene organization between the genomes of these species, a comprehensive pairwise synteny analysis was conducted using MCScan ( Supplementary Figure S2 ). Our analysis revealed extensive synteny as R-genes were largely conserved between species. For example, 88.1% of all the NLR genes in C. inodora (Cin) were syntenic to C. sinensis (Csi), and on average, 53.2% of R-genes in each genome are syntenic with those in another species. On the other hand, the synteny also reflected the accumulation of structural variations, such as rearrangement and translocation, as displayed in the pairs of C. australasica - C. glauca (Cas-Cgl), C. australasica - C. inodora (Cas-Cin) and C. sinensis - C. glauca (Csi-Cgl). It should be noted that no distinct structural divergence was identified between the R-genes in Australian limes and the two cultivated citrus species.

Among all the syntenic R-genes identified in pairwise comparisons, 39 were found to be common across all five species. With these common syntenic R-genes, we constructed a phylogenic tree, which properly grouped the five citrus species into clades according to their origins ( Figure 2A ). In this phylogeny, the two major clades represented the Australian limes and the cultivars. In the Australian lime clade, the two tolerant limes C. australasica and C. indora shared one clade and the resistant lime C. glauca was in a monotypic clade. It’s important to note that this phylogenic grouping was different from the phylogeny based on BUSCO genes, in which C. glauca was positioned in the cultivar clade. In these syntenic R-genes, the numbers of CNL (red), TNL (blue), and NL (green) genes were 10, 10, and 19, respectively ( Figure 2B ). Except for a large number of translocated genes positioned between Chr1 in C. australasica and Chr5 in C. inodora, most of these orthologs were located on the same chromosomes in each species, with nearly all TNL located on Chr3, most NL on Chr5, and most CNL on Chr7. Nearly all the TNL genes were found on Chr3, except one on Chr9 in C. australasica. In contrast, a few NL and CNL genes were located on different chromosomes.

Distinct patterns were found in chromosomal location of R-genes and the synteny of R-genes within each genome. The R-genes are distributed unevenly on the nine chromosomes in the five citrus species, with the majority of NLR genes residing on three chromosomes, Chr3, 5, and 7, in which Chr5 contained the most in each species (75.2-89.6%), followed by Chr3 and Chr7 ( Table 2 ). Up to 75-89.4% of the NLR genes were found in clusters, i.e. at least three genes located within the range of 200 kb. Each of the two cultivars had more clustered NLR genes than the Australian limes. The average cluster size ranged from 7.7 to 10.8 genes, and the largest clusters varied greatly in size across species, ranging from 15 to 74 genes. For the NLR genes not in clusters, they occurred either in pairs or singles, and there were similar proportion of each ( Table 2 ).

On the three R-gene enriched chromosomes, R-genes of the same group tended to form clusters, residing in proximal regions on chromosomes in each species ( Figure 3 ). In particular, nearly all TNL genes were located on Chr3 in tight clusters, whereas most CNL and NLR-genes resided either on Chr5 or 7. In C. clementina and C. sinensis genomes, TNL genes on Chr3 had syntenic genes exclusively on the same chromosome, while TNL genes on the Chr3 of the three Australian species all have syntenic genes on chromosomes other than Chr3. There are more inter-chromosomal syntenic pairs in Australian species compared to the two cultivars, especially C. sinensis, in which only three syntenic NL genes were found between Chr5 and 7. The tandem repeats of the NLR genes on the same chromosomes indicate tandem duplication, while the syntenies of genes from different chromosomes indicate duplication due to transposable elements or more complex rearrangements.

To evaluate the effect of selective forces on the evolution of R-genes in these citrus species, we calculated the rates of nonsynonymous and synonymous substitution ratios (Ka/Ks) on each pair of syntenic R-genes between species using DnaSP6 (Rozas et al., 2017). The comparisons were made in three groups, with the first group (CC) including the two cultivars only, the second group including cultivars versus Australian species (CW), and the third group including the Australian species only (WW) ( Figure 4A ). In each group, the majority of the Ka/Ks ratios were less than 1, indicative of purifying selection. In the CC group, 12 genes have Ka/Ks ≥ 1, with the highest ratio reaching 2.26. In the CW group, there were 25 genes with Ka/Ks ≥ 1, with the comparison between having the most in C. sinensis vs C. inodora. There were seven genes with Ka/Ks ≥ 1 in the WW group.

The Ka/Ks ratios between the two cultivars ranged from 0 to 2.26, with the median being 0.495, and there are 12 genes that have Ka/Ks ≥ 1. In the comparison between cultivars and Australian limes, the median of the Ka/Ks ratios ranged from 0.41 to 0.51. There were 25 genes with Ka/Ks ≥ 1, and the highest ratio was 2.21 between C. sinensis and C. inodora. The comparison between Australian species showed seven genes with Ka/Ks ≥ 1, and the medians ranged from 0.41 to 0.45.

Mapping the Ka/Ks values in the CW group on the chromosomes depicted how selection acted differentially on the R-gene evolution between Australian species and cultivars ( Figure 4B ). The Ka/Ks values tended to be higher in the regions near the telomeres. Some regions on Chr5 of each comparison include low Ka/Ks values, suggesting these regions experienced strong purifying selection. Annotation of the 17 genes with Ka/Ks ≥ 1 between cultivars and Australian limes indicated that all the genes were related to ADP binding. Those from C. clementina were generally associated with coiled-coil domains, and the genes from the C. sinensis were all related to 10,13-epoxy-11-methyl-octadecadienoate biosynthesis ( Supplementary Table S6 ).

Using Orthovenn3 (Sun et al., 2023), the collinearity analysis conducted for orthologous clusters on the R-gene protein sequences identified 57 to 71 orthologous clusters in the five citrus species ( Figure 5A ). About 50 orthologous clusters were shared by all five species, representing 92.8% of the total R-genes ( Figure 5B ). Two unique clusters were identified in C. australasica, comprising four genes; one unique cluster was found in C. glauca and C. sinensis, with two or three genes, respectively. No unique clusters were identified in C. clementina or C. inodora.

Gene family evolution analysis was employed to infer how R-genes diverged among these citrus species during the course of their evolution ( Figure 5C ). According to the dated phylogenetic tree on R-genes, the five species shared a common ancestry approximately 3 million years ago (MYA). The Australian limes underwent the first divergence around 2.75 MYA, in which C. inodora separated from C. australasica and C. glauca, which underwent another divergence of about 2.0 MYA. In contrast, the divergence between the progenitors of the two cultivars occurred 1.5 MYA.

All the R-genes in the three Australian limes belonged to one gene family and experienced gene family expansion approximately 3.0 MYA, followed by another gene family change of about 2.0 MYA, including contraction of 6 gene families and expansion of 2 gene families. In contrast, the first gene family change in the two citrus cultivars occurred around 1.5 MYA, in which all 11 gene families experienced contraction.

To examine the composition and distribution of R-gene domains, we retrieved domain sequences and their chromosomal coordinates from TNL genes of each species. A similar composition of R-gene domains was found to reside on these genes across all the species ( Figure 6A ). In most genes, there was one single copy of full-length TIR and varying abundance of LRR domains. A total of six LRR domains were identified, with LRR1, 3, 4, and 8 common to all five species while LRR5, and LRR6 occurred rarely. The average numbers of the LRR4 domain (Pfam ID: PF12799) per gene were higher in the three Australian limes (2.9-3.2) compared to each of the two cultivars (2.5-2.6). The number of LRR8 (Pfam ID: PF13855) was lowest in C. clementina, and similar in the other four species. Domains of LRR1 and LRR3 only occurred in single copies on a gene, and there were, on average 0.29 and 0.83 per gene across the species.

To facilitate visualization and comparison of the chromosomal distribution of R-gene domains, we reduced the gene number in each species down to 20 using Treemer (Menardo et al., 2018), which evaluates the redundancy of phylogeny and only keeps the nodes that contribute the most to the phylogenetic diversity ( Figures 6B–F ). While TIR domains were generally located near the N-terminals in each gene, the LRR domains are 400-600 bp away, with most occurring in tandem arrangement. LRR4 and LRR8 are the dominant domains compared to the other two, and their copy numbers varied highly between genes, ranging from one to as many as 10. The distribution and organization of the LRR domains were highly variable in each species, and no clear patterns were demonstrated between resistant and susceptible types.

The LRR domains in R-genes are often involved in specific recognition of pathogen effector molecules. Among the eight LRR domains identified in the five citrus species, four LRR domains (LRR1, 3, 4, and 8) were found to be common in each species. To facilitate comparison of the LRR domains between species, we employed EMBOSS-CONS (Madeira et al., 2024) to generate consensus protein sequences of each domain in each species ( Figure 7 ). The alignment of the consensus sequences revealed that LRR1 and LRR3 are relatively short and highly conserved across the five species, with most positions residing with identical residues ( Figures 7A, B ). In contrast, the sequences in LRR4 are highly divergent, and there are many positions with non-consensus residues or gaps ( Figure 7D ). Except for C. sinensis, the other four species harbor several long stretches (≥5) of insertions composed of non-consensus residues missing in some species. In addition, the sequences of conserved positions are short and frequently interrupted by non-consensus residues. In the three Australian limes, the first positions are conserved, which all have an N-residue, but a gap and a serine residue in C. clementina and C. sinensis, respectively ( Figure 7D , highlighted). The LRR8 domains are moderately conserved ( Figure 7C ). The two cultivars have two consensus residues at positions 8 and 11, which are non-consensus in the three Australian limes ( Figure 7C , highlighted).

Motif analysis detects important molecular features, such as nucleotide binding sites, and predicts protein interaction domains on the target sequences. In this study, we identified the top five motifs from each R-gene, and most of these motifs are about 20 amino acids long ( Figure 7E ). The alignment of these motifs indicates the first domain was highly conserved among all the species, and the second motif is conserved in large regions of all five species. In contrast, no conserved segments were found in the other three motifs. The first two motifs are related to NBARC function, while the other three were not associated with any known function. According to concatenated motif sequences, phylogeny indicated that C. sinensis was closer to the Australian species than C. clementina.

In this study, we used a recently developed pipeline FindPlantNLR to comprehensively identify and annotate NBARC type R-genes from the five citrus species, including three wild Australian limes and two cultivated citrus species. This pipeline produces highly robust and reliable R-gene identification and annotation, enabling R-gene retrieval in unexplored genomic regions (Chen et al., 2023). As demonstrated in this study, FindPlantNLR uncovered many R-genes in each species that were unidentified in the original genome annotation. The proportion of the newly identified R-genes reached as high as over 70% in C. australasica and C. glauca and over 30% in other three species. In general, the NBARC loci cannot be accurately predicted by the automated gene annotation pipelines due to the repeat masking, a necessary step during genome assembly and gene annotation to avoid local genome assembly collapse and annotation errors (Tørresen et al., 2019). Since NBARC genes are typically organized in clusters of tandemly duplicated sequences, the R-gene loci are often obscured during genome annotation and inadvertently excluded from detection with the use of R-gene search pipelines that solely rely on predicted genes derived from genome annotation (Meyers et al., 2003; Andolfo et al., 2013). FindPlantNLRs overcomes this limitation by directly annotating the genome and may explain the discovery of high numbers of previously masked R-genes in this study. Our results show that this or similar pipelines that leverage direct annotation in unmasked genomes is essential for accurate and comprehensive prediction of the R-genes in plants.

In the two cultivars C. sinensis and C. clementina, we identified 689 and 761 NBS genes from the genomes, respectively, numbers higher than the previously reported R-genes. For example, only 111 NBARC domains were identified in protein sequences of C. sinensis (Yin et al., 2023). Combining open reading frames (ORFs) search and protein sequences selection resulted in prediction of higher numbers of NBARC genes in C. sinensis (Wang et al., 2016) than the earlier report (Yin et al., 2023), but still significantly lower than our prediction. Further, only about 50% of NBARC genes in the study by Wang et al. (Wang et al., 2015) could be classified as NLR genes, remarkably lower than our classification, ranging at 84-85%. Protein sequence comparison using PLASTP revealed that each NBARC gene identified in these studies was within the complement of R-genes predicted in our study in both species ( Supplementary Table S3 - S5 ). Considering the use of automatically annotated genomes in these studies, the above results reinforced the notion that approaches of direct genome detection are essential for more comprehensive and complete annotation of R-genes.

Comparing the overall R gene complements, the two cultivars harbored more NLR genes than the Australian limes, especially C. inodora and C. glauca ( Table 1 ). It is not uncommon that the abundance of NBS genes in cultivated species is higher than in their wild relatives. For example, the numbers of NLR genes in the Asian cultivated rice (Oryza sativa L.) were found to be substantially higher than those in the wild ancestors (Mizuno et al., 2020), and cultivated soybean showed roughly 3-fold more TNL than its wild relative (G. latifolia) (Liu et al., 2018). Such R-gene expansion and duplication may have resulted from domestication and cultivation. Similar effects might also be at play in cultivated citrus, which are either natural hybrids or are bred mainly through interspecific hybrid and/or admixture (Wu et al., 2018). It has been suggested that retention of duplicated R-genes often incurs fitness costs, and as such natural selection tends to maintain a limited number of resistance genes in favor of more significant growth and reproduction. At the same time, domestication often enriches R-genes (Barabaschi et al., 2020).

Nevertheless, it should be noted that the abundance of R-genes in the genomes may not be directly related to resistance capability. As observed in this study, the Australian limes were found to have fewer NLR-genes than the cultivated species, though they have proved to be more tolerant/resistant to HLB (Ramadugu et al., 2016b). In the R-gene mediated defense system, the constitutive expression of a core set of R-genes is essential in implementing on-going defense status (von Dahlen et al., 2023), and different genotypes may vary extensively in the basal expression of most R-genes. For example, in a survey of 45 gene families studied in 19 accessions of A. thaliana, two R-gene subfamilies were found to be among the top three families of highly expressed genes, and the extent of differential expression for R-genes was surprisingly high, reaching up to 350-fold difference between accessions (Gan et al., 2011). Though R-genes are believed to function in a gene-for-gene manner, their expression patterns are also shaped by evolutionary paths and subjected to the influence of environmental factors (MacQueen and Bergelson, 2016), which may also explain the lack of association between R-gene abundance and resistance levels.

One striking feature of plant NBS-LRR genes are their frequent clustering arrangements. Our analysis showed NLR genes in each citrus species are preferentially mapped to three chromosomes (3, 5, and 7) and mainly occurred in clusters ( Table 2 and Figure 3 ). This result is consistent with previous findings in C. sinensis and C. clementina, in which 76.9 and 84.9% of the respective NBS genes were found in clusters (Wang et al., 2016). The high percentages of clustered NBS genes were also reported in other species such as rice (Yang et al., 2006), Arabidopsis (Guo et al., 2011), grapevine, and poplar (Yang et al., 2008). The synteny analysis indicated that the clustering of R genes may have arisen from both tandem duplication and translocation, which resulted in clustered sequences along or across chromosomes, respectively ( Figure 3 ). The clustering arrangement is believed to be advantageous in providing coregulatory benefits and a broader detection spectrum against pathogens (van Wersch and Li, 2019). Accumulating evidence suggests the NLR genes continue to evolve in complexity and tend to function in higher-order configurations, such as NLR pairs or networks, formed by clustered NLRs, rather than dispersed genes or singletons (Contreras et al., 2023). For example, it has been found that coupled NLR genes are required to initiate resistance against a single pathogen in several species, such as genes associated with viral resistance in Arabidopsis thaliana and tobacco, leaf rust resistance in wheat, and blast resistance in rice (Zhai et al., 2014). The genetically clustered NLR genes in tight physical proximity supply raw genetic material for the acquisition of new resistance in the processes of sub-functionalization or neo-functionalization (Michelmore et al., 2013). The Australian wild limes included in the study originated from Australia where HLB disease or the associated pathogens have not been reported.

Phylogenetic relationships of R-genes depict the evolutionary connections of disease resistance between species. Using the 39 NLR genes that were common to each citrus species, we demonstrated that the phylogenetic inference precisely reflected their HLB-resistance levels, i.e. the three clades corresponded to the three resistance categories, i.e. resistant (C2) (containing C. glauca), tolerant (C3) (containing C. australasica, and C. inodora), and susceptible (C7) (containing C. clementina, and C. sinensis), based on Ramadugu’s (Ramadugu et al., 2016b) category system ( Figure 2 ). This R-gene based phylogenetic inference was more accurate than the one using BUSCO genes ( Figure 1B ), which generated low support values on each branch and wrongly grouped an Australian lime (C. glauca) with the cultivated citrus species, suggesting R-genes may have evolved at a faster rate than the BUSCO genes. However, our domain structure and motif analyses ( Figures 6 , 7 ) revealed no distinct patterns unique to either group, indicating resistance-related genomic characters may lie beyond the domains or motifs; the regulatory components may play an important role in determining resistance.

It has been demonstrated that R-gene expression is controlled by a complex yet tight regulatory system (Stokes et al., 2002; Holt et al., 2005; Li et al., 2007; Huot et al., 2014). The regulation of R-gene expression can operate on multiple levels, including transcriptional and epigenetic regulation, RNA interference (RNAi), splicing and translation, and post-translational controls (Kapos et al., 2019). In addition to the expression in the presence of pathogens, some R-genes were expressed variably across species irrespective of infection status (von Dahlen et al., 2023). Also, R-gene expression is highly tissue-specific, as evidenced in several transcriptomic studies (Chen et al., 2007; Zhai et al., 2014; Sharma et al., 2017). In addition, R-genes may be activated in response to environmental factors alone without biotic stimuli (MacQueen and Bergelson, 2016). Given the complexity of R-gene regulation, the HLB resistance observed in Australian limes may mainly be derived from regulatory mechanisms (in addition to the presence of essential R-gene domains).

The Ka/Ks ratios, which measure the relative impact of diversifying and purifying selection on R-gene sites, were used to estimate which sites in the R-genes from the Australian limes were advantageous over cultivated citrus species. In the five Australian limes, we identified a total of 25 R-genes that had Ka/Ks ratios greater than one, indicative of positive selection in these sites. The low number of sites with high Ka/Ks ratios indicates limited evolutionary pressures in the Australian limes, suggesting the high conservation of the orthologous gene pairs between the two groups. Further, among these R-genes, 20 were located in clusters, and five existed as singletons. This finding is consistent with the proposal that R-gene singletons tend to maintain sequence stability and functional conservation under strong purifying selection, whereas R-genes in clusters undergo fast evolution to facilitate functional innovation in coping with changing pathogenic threats (Zhang et al., 2019). Similarly, in comparing cultivated citrus species, the median Ka/Ks rate ratio of singletons is less than those in clusters (Wang et al., 2015). Together, the R-genes identified in Australian limes may have accumulated beneficial mutations and may potentially serve as molecular markers to assist in breeding for pathogen-resistant cultivars.

The polymorphisms present within R-genes are subjected to the evolutionary forces. Fossil and molecular evidence has suggested that the genus citrus originated in southeast Asia, approximately 8 million years ago (Xie et al., 2023), and Australian citrus species arose during a major ancient species dispersal (Australian radiation) approximately 4 MYA (Schwartz et al., 2015; Wu et al., 2018). Chloroplast genome phylogeny (Wu et al., 2018) and pangenome analyses (Huang et al., 2023) indicated C. glauca diverged from C. australasica between 2-4 MYA. Our orthologous and phylogenetic analysis using NRL genes reflected the divergences between these species ( Figure 5C ), but on a relatively smaller time scale. This discrepancy may be due to the conservation in the R-genes, which are unlikely to show lineage sorting during rapid radiation and speciation events. The high overlapping of R-gene clusters also confirmed conservation among the R-genes in these species ( Figure 5A ). It is noticeable our result indicated the divergence between the progenitors of two cultivars were inferred to be 1.5 MYA, though selection of C. clementina was a recent event. This reflects some R-genes may have evolved at a faster rate, and thus inflated the estimation of evolutionary timeline between these two species.

The C-terminal LRR domains mediate pathogen recognition in NLR proteins through protein-protein interactions. As such, LRR domains contribute the most to R-gene polymorphism, as evidenced in the analysis of sequence mutations, protein secondary structure, and three-dimensional structures (Ratnaparkhe et al., 2011). In this study, we identified six LRR domains in the five citrus species, which exhibited high variation in abundance and genomic distribution ( Figure 6 ). This result is consistent with the previous study in three citrus species, in which the LRR motifs showed high variation in sequences and repeat numbers (Wang et al., 2015). The sequence diversity in the LRR domains is consistent with the role of LRR domains in constituting a sensor domain that interacts with various molecular partners in detecting a variety of ever-evolving pathogens (Takkouche et al., 2023). Noticeably, LRR4 domain appeared to be more polymorphic between species compared to the other three LRR domains ( Figure 7 ) and was more abundant in the three Australian limes than in each of the two cultivars ( Figure 6A ). Usually, one LRR domain detects a specific target effector, but plant LRR proteins must undergo conformational changes to induce the downstream defense responses (Liu et al., 2023). Therefore, the relative abundance of LRR4 in the Australian species may allow for expanded conformational complexity in mediating pathogen recognition, thus increasing regulatory capacity in coping with pathogenic attacks.

In this study, we analyzed the genomic complement of R-genes in five citrus species to characterize the differences between Australian limes and cultivated species, thereby paving the way for the development of tools for genome-assisted breeding for HLB-resistant varieties. The syntenic analysis indicated the R-genes sequences contributed to the difference in HLB-resistance levels. However, substantial similarities in the genomic structure of R-genes were revealed in the five citrus species, and the identified polymorphisms were insufficient to distinguish between the two groups. These findings suggest that the HLB resistance in Australian limes may involve mechanisms other than R-genes. As suggested in a transcriptomic study (Weber et al., 2022), the resistance mechanisms in C. australasica may include phloem callose formation, redox control, phytohormone mediated signaling, secondary metabolites, secretion of pathogenesis-related (PR) proteins (i.e., Cys-rich secretory proteins and PR1-like proteins). Future investigations depicting association of HLB resistance and hybrids between resistant and susceptible citrus species in combination with RNA-seq data would yield more insights into R-genes along with other mechanisms that are responsible for HLB-resistance.