Introduction Resistance (R) genes confer to plants innate immunity against a broad spectrum of pathogens, including viruses, bacteria, fungi and nematodes, through expressing matching avirulence genes in a “gene-for-gene” manner 1 – 5 . There are five main recognized groups of R proteins based on differences in the domains of cloned R genes 1 , 6 – 10 . Among these, the largest group encodes receptor-like proteins that feature a conserved nucleotide-binding site (NBS) domain accompanied by a C-terminal leucine-rich repeat (LRR) domain that can recognize antigen proteins of a pathogen. The NBS R genes are further separated by an N-terminal Toll/Interleukin-1 Receptor (TIR domain) or a coiled-coil (CC motif) homology region 11 – 13 . The TIR domain is the component of a large R gene subgroup in dicot, conifer, and moss genomes, but the domain is almost nonexistent in monocots 8 , 12 , 14 , 15 . Whole-genome duplication (WGD) and tandem duplication (TD) contributed to the expansion of the largest NBS R gene group in plant genomes 16 . The NBS R genes have been studied via whole-genome sequencing in dicot and monocot species, including Arabidopsis thaliana, Oryza sativa, Populus tremula, Carica papaya, Medicago truncatula and Vitis vinifera 7 – 9 , 17 – 19 . According to previous studies, approximately 0.2–2% of the genes in plant genomes are NBS R genes, with papaya—which has not experienced recent whole genome duplication—retaining relatively fewer NBS R genes in its genome 17 . In the genomes of A. thaliana, O. sativa, and C. papaya, the majority of the NBS R genes are clustered in local regions, which is consistent with the consequences of tandem duplication events 7 – 9 , 16 , 17 , 20 – 22 . R genes play an important role in monocots, including field crop species. Previous studies have mapped and cloned many R genes involved in disease resistance and have found that some R genes were under strong selection related to pathogen Avr genes 23 . This work has contributed to the breeding of disease-resistant crops in agricultural production. Among the more well-studied R genes, the majority are NBS R genes. In Zea mays, Rxo1 has been cloned as a resistance gene against bacterial streak disease 24 . Rxo1 is an NBS R gene that confers resistance to a set of pathogen strains, including an unrelated pathogen that causes bacterial stripe in Sorghum bicolor and Z. mays. Rxo1 retained its resistance function when used as a transgene into O. sativa as a defense against the bacterial streak disease, demonstrating the shared functional mechanism and the feasibility of cross-species utility of R genes among grass species. In O. sativa, more than 20 NBS R genes have been cloned to confer resistance against the blast and bacterial blight diseases. Fei et al. 25 , Zeng et al. 26 , and Wu et al. 27 independently cloned the blast resistance gene Pia from the loci of NBS R gene TD arrays by genetic mapping of resistant cultivars of O. sativa. Guo et al. cloned several O. sativa blast resistance genes through paralogous comparison based on sequence homology 28 . Furthermore, Yang and co-authors analyzed rapidly evolving R gene families in maize, sorghum, brachypodium, and rice, and found that after transforming some of these genes into rice, the transformed plants gained blast disease resistance 29 . This suggested that fast evolving R genes play an important role in functional innovations of resistance. Though many NBS R genes have been determined to maintain resistant traits, there are still a large number that have not been functionally characterized and thus are a rich resource for resistance function innovation, mechanism investigation, and potential breeding utilization. Systematic determination and comparison of R genes among multiple closely related species will help to reveal evolutionary features with regard to expansion, divergence, and function innovation. R genes have been studied in many grass species, including O. sativa, Z. mays, and Setaria italica 9 , 13 , 22 , 30 following whole-genome sequencing of these species, and several studies have been performed with limited comparisons among 3–4 species 31 , 32 . However, genome-wide determination and comparison of R genes across all sequenced grass species has not yet been done. In this study, we focused on all grass species whose genomes have been sequenced and released (11 genomes), together with two outgroup species. We performed comprehensive comparisons of R genes after genome-wide determination under the same rules in these 13 species. We found that R gene loci showed high levels of conservation among species, and R genes from TD and singleton loci adopted different strategies in the functional evolution of disease resistance. Results Identification of R genes in Setaria italica and closely related grass species Genome sequences and predicted gene information for 13 species (Table 1 ) were downloaded from the databases phytozome and PlantGDB for R genes analysis. To check their phylogeny relationships accurately—as a guide for analyzing R gene evolution—we performed pairwise genomic synteny analysis between S. italica and each of the other species (Methods). Using the genome of S. italica as a control, P. virgatum and Z. mays were confirmed as tetraploids; T. aestivum is hexaploid, while the other species are diploids (the distant outgroup species M. acuminate and A. thaliana are not applicable) (Supplementary Figs S1 – S9 ). Based on the syntenic gene datasets, we obtained 357 groups of complete syntenic orthologs (syntenic orthologs exist in all studies genomes) among the 12 grass species (except for the dicot A. thaliana). We further extracted 10,443 K s loci (coding variants that would not cause amino acids changes to the corresponding proteins) from these 357 syntenic ortholog groups. These loci were then concatenated and used to build the phylogenetic tree under the neighbor-joining algorithm in MEGA. As shown in Fig. 1 , with M. acuminate as the outgroup, the evolutionary relationships of these 12 grass species were clear, meaning that the phylogeny serves as a good reference for evolutionary analysis of R genes among these species. Table 1 Genomic information of 12 grass species and one outgroup species studied in this work. Index Species #Chr* #Genes Genome (Mb) WGP** 1 Setaria italica 9 34584 405.74 1 2 Setaria viridis 9 35214 394.9 1 3 Panicum virgatum 18 91838 1689.57 2 4 Panicum hallii 9 37232 554.13 1 5 Sorghum bicolor 10 34211 732.15 1 6 Zea mays 10 88760 2067.86 2 7 Oropetium thomaeum / 28437 243.17 1 8 Triticum aestivum 21 107891 14547.26 3 9 Brachypodium stacei 10 29898 234.14 1 10 Brachypodium distachyon 5 34310 271.16 1 11 Oryza sativa 12 42006 374.47 1 12 Musa acuminate 11 36528 472.96 — 13 Arabidopsis thaliana 5 27751 119.67 — * Chr: Chromosomes; **WGP: whole genome polyploidization. Figure 1 Phylogenetic tree built on K s loci of syntenic orthologous genes among the 12 grass species. Red stars denote whole-genome duplication events; double red stars denotes whole genome triplication events. Domain compositions of well-determined R genes in model species A. thaliana and O. sativa were used to search for R genes in those 13 genomes. A total of 202 and 456 R genes of A. thaliana and O. sativa, respectively, were collected from previous studies 8 , 22 . First PfamScan 33 was used to search for domains in these R genes. In all, 55 domains were identified from R genes of the two model plant species (Supplementary Table S1 ). These domains were further used to determine candidate R genes in the genomes of all 13 species using the tool hmmsearch 34 . After that, we obtained 1,264 R gene candidates in S. italica. By filtering genes without the NBS domain, we finally obtained 535 NBS R genes (referred as R genes hereafter) in S. italica. Figure 2 shows the phylogenetic relationships of these 535 R genes determined in S. italica using the maximum likelihood algorithm in MEGA. With the same method, we determined a set of R genes ranged from 97–2,747 copies (accounting for 0.34–2.55% of the total number of genes) in these 13 genomes (Table 2 ). The hexaploid T. aestivum had the most R genes (2,747), while the tetraploid P. virgatum had the second most R genes (1,267). In the diploids, O. sativa (587) and S. italica (535) had the most R gene copies, while O. thomaeum had the least number of copies (97). The R gene copy numbers show association to the polyploidization levels of the corresponding species, except for the tetraploid Z. mays, which has relatively fewer (306) R genes. When comparing numbers of R genes to total genes, T. aestivum had the highest ratio of R genes (2.55%), while O. thomaeum had the smallest ratio (0.34%), followed by Z. mays (0.35%), M. acuminate (0.52%), and A. thaliana (0.83%). R genes can be further classified into three main groups as NBS, CC-NBS, and TIR-NBS. The 535 R genes in S. italica can be further separated into 14 subgroups on the basis of composition and order of domains (Table 3 ). Among these R genes, the NBS group R genes had the most copies, while the TIR-NBS group R genes had the least number of copies. The other genomes shared similar distribution patterns of R gene subgroups, except for A. thaliana, in which TIR domain-containing R genes had the most copies (Supplementary Table S2 ). In contrast to A. thaliana, the 12 grass species had very limited numbers of copies of TIR-containing R genes, which is consistent with previous reports 8 , 12 , 14 , 15 . Figure 2 The phylogenetic tree of 535 R genes in S. italica. The tree was built with the maximum-likelihood algorithm by MEGA software. The red full circles denote S. italica specific R genes, i.e., S. italica R genes that have no syntenic orthologs in any of the other 11 grass species. Branches colored in orange are R genes from tandem duplication arrays, others are R gene singletons. Table 2 Statistics summary of R genes in the 13 studied species. Index Species #R candidate #NBS R NBS R/Gene #R TD* #TD R** TD R/NBS #R singleton #R locus 1 S. italica 1,264 535 1.55% 106 303 0.57 232 338 2 S. viridis 1,152 465 1.32% 92 241 0.52 224 316 3 P. virgatum+ 2,949 1,267 1.38% 212 524 0.41 743 955 4 P. hallii 1,101 420 1.13% 86 234 0.56 186 272 5 S. bicolor 1,036 425 1.24% 82 233 0.55 192 274 6 Z. mays+ 1,101 306 0.34% 36 79 0.26 227 263 7 O. thomaeum 410 97 0.34% 4 8 0.08 89 93 8 T. aestivum++ 5,712 2,747 2.55% 547 1,588 0.58 1,159 1,706 9 B. stacei 945 377 1.26% 78 203 0.54 174 252 10 B. distachyon 1,091 464 1.35% 105 267 0.58 197 302 11 O. sativa 1,390 587 1.40% 116 339 0.58 248 364 12 M. acuminate 863 190 0.52% 24 67 0.35 123 147 13 A. thaliana 896 229 0.83% 49 124 0.54 105 154 *R TD means tandem duplication array that is composed of R genes; **TD R means the R genes comes from a tandem array; +: tetrapolyploidization; ++: hexapolyploidization. Table 3 The subgroups of 535 R genes in S. italica. Index Domains #Copy Sum 1 CC-NBS-LRR 159 230 2 CC-NBS 41 3 CC-NBS-NBS-LRR 14 4 CC-NBS-NBS 12 5 CC-NBS-LRR-LRR 3 6 CC-LRR-NBS-LRR 1 7 NBS-LRR 156 302 8 NBS 88 9 NBS-NBS 45 10 NBS-NBS-LRR 10 11 NBS-LRR-LRR 2 12 LRR-NBS-NBS-LRR 1 13 TIR-NBS 2 3 14 TIR-NBS-TIR 1 Total 535 High level of synteny relationships of R loci among grass species There are two groups of R genes based on whether or not they are distributed continuously or occur singly. The closely distributed R genes are called tandem R genes (R genes from TD arrays), while the others are called R gene singletons. In order to separate the two groups of R genes, we considered continuously distributed R genes that were interrupted by less than five non-R genes as R genes in TDs (referred as R TDs hereafter) (Methods). In all, we determined 106 R TDs (comprising 303 R genes) in S. italica, accounting for 56.63% of a total of 535 NBS R genes (Table 2 ). The remaining 232 R genes are singletons. These R TD genes are colored in orange in Fig. 2 . As shown in Fig. 2 , the R TD and singleton genes are generally mixed within the phylogeny, with more R singletons clustered in two more divergent branches (bottom-right part of Fig. 2 ). This indicates that these R singletons were formed much earlier and were conserved as singletons for a longer time than R TDs in the upper branches of the tree. Using the same method, we obtained 4–547 R TDs (comprised of 8–1,588 R genes) in the other genomes, accounting for 3.72–57.75% of the total number of R genes in corresponding genomes (Table 2 ). We further selected one of the R genes in each of the 106 TDs in S. italica as a representative gene of each R TD. Together with the 232 R singletons, there were 338 R gene loci (referred as R loci hereafter) in the genome of S. italica. Similarly, we determined 93–1,706 R loci in the other genomes (Table 2 ). Additionally, we investigated the domain-related subgroups of R genes in R TDs, and found that different R domain subgroups showed no significant distribution differences in relation to R TDs. However, for the subgroups of R genes that had less than five copies, we found that they had a high probability to be from R TDs. These R genes with limited copies show duplicated or abnormal orders of domains. It seems reasonable that those genes are R genes from TDs that originated from domain rearrangement of duplicated R genes in the same TDs. Most R genes in S. italica had syntenic orthologs in the other grass genomes. Based on whole-genome synteny analysis (see Methods), we obtained syntenic orthologous gene sets between each pair of the 13 species, i.e., 78 pairwise comparisons in total (Supplementary Table S3 ). Based on the datasets, we checked the syntenic orthologs of the 338 S. italica R loci in other genomes. We found that 315 out of 338 R loci (93.20%) had syntenic orthologous loci in at least one of the other genomes, and 131 R loci (38.76%) had syntenic orthologs in more than half (6) of all compared genomes (Supplementary Fig. S10 ). We further analyzed the synteny relationships of R loci for the two R gene groups separately. For these 106 R TDs in S. italica, we found that 102 (96.23%) had syntenic orthologous loci in at least one of the other species, among which 94 loci were R TDs in at least one of the other species; for the 232 R singleton loci, we found that 213 (91.81%) had syntenic orthologs in other species, among which 98 syntenic loci were R TDs in at least one species. These four R TDs and 19 R singletons that are specific to S. italica are denoted by red circles in Fig. 2 . There were no clear clustering differences to other R genes in S. italica, indicating that they originated from transposition duplications of these R genes that have synteny relationships to other genomes. R genes evolved rapidly through tandem duplications and are conserved in singletons Most of the R loci retain syntenic relationships to R loci in the other grass genomes. After removing the two distant outgroups M. acuminate and A. thaliana, which have limited syntenic orthologs to the 11 grass genomes, as well as T. asteria whose genome was extremely expanded, we obtained syntenic orthologous gene families among the ten grass species based on the aforementioned pairwise syntenic gene datasets. In all, we obtained 112,723 such syntenic gene families among the ten genomes. We found that 785 out of the 112,723 syntenic orthologous families belonged to the R loci syntenic gene families. Among these, 359 syntenic families of R loci had syntenic orthologs in two or more species, while the other 426 R loci families were unique in only one of the ten grass genomes (17–143 R loci families were species-specific to one of the 10 genomes) (Supplementary Table S4 ). For the 359 syntenic families, 168 families contained R TDs in at least one of these syntenic R loci (46.90%). For the 426 species-specific R loci, 85 were R TDs (19.95%), which is significantly lower than that of shared R loci among species (χ 2 test, P = 1.63 × 10 –8 ). This indicates that the shared R loci tended to retain more R TD sequences in the evolutionary expansion of R genes, while the species-specific R genes tended to be retained in fewer copies as R singletons. R genes in TDs apparently evolve more rapidly than R singletons. We calculated K a (the rate of non-synonymous mutations per non-synonymous locus) and K s (the rate of synonymous mutations per synonymous locus) values between pairs of genes from the R syntenic orthologous families among the ten grass species in order to estimate the selection strength imposed on syntenic R genes. We further investigated the frequency distributions of these K a , K s and K a /K s values for R genes from TDs and singletons separately. For R TDs, we have 39,855 such pairwise syntenic orthologs comparisons, while for R singletons, we obtained 1,598 pairwise comparisons. The K a values from R gene singletons were much less than those of R genes from TDs (Supplementary Fig. S11 ). This suggests that R genes in TDs have been evolving and diverging much faster than R gene singletons. We further examined the differences of K a /K s values within the same comparisons; the analysis showed that R gene singletons have smaller K a /K s values than those of R genes from TDs (Fig. 3 ). The result indicates that the R gene singletons are under much stronger purifying selection than R TDs genes. These findings together suggest that R genes have expanded and rapidly evolved through TDs, while R genes with important original or newly evolved functions would be conserved as or reduced into singletons. Figure 3 The frequency distribution of K a /K s values between syntenic R gene pairs from R singletons (colored in red) or R TD arrays (colored in light orange). R genes are under selection in populations of S. italica R genes have accumulated more mutations than other genes. From the whole-genome resequencing data of 916 S. italica accessions reported and released previously 35 , we determined 2,156,941 SNPs and 189,440 InDels variants from the S. italica population after mapping the resequencing data to the reference genome of S. italica 36 (See Methods). We then annotated these variants (detailed information is listed in Supplementary Tables S5 and S6 ). There were 704,436 SNPs and 89,804 InDels located in the genic regions of 34,584 genes in S. italica. We further considered SNPs that caused non-synonymous mutations or splicing alterations (within two bps of splicing loci) as well as InDels located at coding sequences or splicing loci as functional variants, as they lead to changes in the final protein sequences. There were 79,296 such functional variants obtained (70,671 SNPs and 8,625 InDels), an average of 2.29 variants per gene. We then examined the 535 R genes and found that they had accumulated 9,006 functional mutations (8,466 SNPs and 540 InDels), which is approximately16.83 functional variants per R gene. This was significantly higher than the average level of all genes (2.29 functional variants per gene) (χ 2 test, P &lt; 2.20 × 10 −16 ) (Table 4 ). Furthermore, we checked the two aforementioned R gene groups and found that the R genes from R TDs accumulated 20.70 functional mutations per gene, while these singletons only had 11.78 such variants per gene, significantly lower than that of R TDs (χ 2 test, P = 3.68 × 10 −10 ). This suggests that R genes in tandem arrays have evolved more rapidly than R gene singletons. Table 4 Functional mutation comparison between total genes and R genes. #Genes #Functional mutations Functional mutations/Genes P values (χ2 test) Total gene 34,584 79,296 2.29 &lt;2.20 × 10 −16 NBS R gene 535 9,006 16.83 R singleton 232 2,734 11.78 3.68 × 10 −10 R TD 303 6,272 20.70 R genes are under stronger selection than other genes. Using the variome datasets from S. italica populations, we evaluated the selection pressure imposed on genes by calculating the measures π 37 and Tajima’s D 38 . The values of π range from 0 to 1 and estimate the level of polymorphism of given sequences in a population, while Tajima’s D estimates the degree of purifying or balancing selection. We obtained 34,584 values of π and Tajima’s D corresponding to 34,584 genes in the genome of S. italica. Figs 4 and S12 show the distributions of these values for all genes and R genes, respectively. The R genes are distributed as a steep peak around the π value of 0.35, which suggests that they are under more divergent selection compared to all genes in the population of S. italica. As mentioned above, we estimated the selection strength for groups of R TDs and R singletons separately and found that R singletons exhibited distribution peaks at lesser values of π and Tajima’s D (Figs 4 and S12 ). This suggests that the R genes from TDs are under balancing selection, while areas of R singleton genes are under purifying selection. Taken together, the results suggest that R genes have undergone tandem duplications that promoted rapid evolution, while singletons enabled R genes to retain and conserve sequence and function. Figure 4 The distribution of π values of total genes (colored in black) and R genes (light blue) in the population of S. italica. R genes were further separated into R singletons (red) and R tandem arrays (light orange). Newly mapped R genes show strong divergence among species Syntenic orthologs of newly mapped functional R genes are under divergent selection among grass species. Pi37 is an R gene that confers the trait of blast disease resistance in rice 25 . We examined the synteny orthologs of Pi37 in the 11 grass genomes and found that ten of the 11 species had syntenic R gene loci to Pi37 except Z. mays (Table 5 ), with two loci in S. italica and O. thomaeum being involved in micro-translocations in local genomic regions. Although these syntenic orthologs were inherited from a common locus of a most recent common ancestral genome, their coding sequences have heavily diverged (identity: 27.54–68.34%, coverage: 38.84–100%; Table 5 ). This indicates that these syntenic R loci have rapidly evolved and diversified in different species, and that they may have evolved different resistance functions. The previously reported blast resistance gene Pi37 is located in an R TD (syntenic loci in P. hallii and T. aestivum also are R TDs) in O. sativa. This R TD locus may have originated after the divergence of these species, and the R TD then contributed to blast resistance functional innovation in O. sativa. The syntenic gene of Pi37 in S. italica is Seita.5G337400.1, which is an R singleton and is involved in a local translocation. Seita.5G337400.1 showed the lowest identity to Pi37 of O. sativa, and thus it may not have the same blast resistance function as that of O. sativa Pi37. Meanwhile, Seita.5G337400.1 has been under selection in populations of S. italica as indicated by π (0.18) and Tajima’s D (0.20) compared to other genes (Figs 4 and S12 ). However, the target of selection on Seita.5G337400.1 in S. italic populations should be different to that of blast resistance of Pi37 in O. sativa. Table 5 The syntenic orthologous genes of O. sativa blast disease resistance gene Pi37 in the other 11 grass species. Species Pi37 synteny orthologs Identity% Coverage% E-value Tandem array O. sativa Os01g57310.1 (Pi37) / / / Os01g57270.1; Os01g57280.1; Os01g57310.1; Os01g57340.1 B. distachyon Bradi2g51807.3 68.34 100 0 / B. stacei Brast01G124500.1 68.03 100 0 / T. aestivum TraesCS3A02G316800.1 66.26 56.90 0 TraesCS3A02G316800.1; TraesCS3A02G316900.1 TraesCS3B02G349500.1 66.29 47.60 0 TraesCS3B02G349500.1; TraesCS3B02G349600.2; TraesCS3B02G349800.1 TraesCS3D02G314300.4 64.56 55.12 0 TraesCS3D02G314300.4; TraesCS3D02G314400.1 O. thomaeum Oropetium_20150105_12744A 45.11 38.84 4.42E-84 / Z. mays / / / / / S. bicolor Sobic.003G317300.2 62.11 100 0 / P. hallii Pahal.E01385.1 61.95 100 0 Pahal.E01385.1; Pahal.E01390.1 P. virgatum Pavir.5KG572900.1 61.53 92.48 0 / S. viridis Sevir.5G344900.1 55.43 100 0 / S. italica Seita.5G337400.1 27.54 48.14 1.68E-40 / We further analyzed the syntenic orthologous genes of Z. mays Rxo1 in the 11 grass species. Rxo1 is an R gene that has been mapped and reported to confer resistance against the bacterial streak disease 24 . There were seven species that retained such syntenic loci (Supplementary Table S7 ), and all are involved in local micro-translocations. We found that all of these syntenic loci were R singletons, except for Z. mays Rxo1 located in an R TD. The loci show even poorer sequence similarity (identity: 26.55–27.70%, coverage: 38.84–100%) than the Pi37 syntenic orthologs, indicating greater divergence among these syntenic orthologous genes of Rxo1. The syntenic gene of Rxo1 in S. italica is Seita.6G072000.1, which has been under diversifying selection in S. italica (π: 0.38 and Tajima’s D: 1.48). These findings for Pi37 and Rox1 supported the suggestion that R TDs contributed to the functional innovation of R genes in different species to cope with distinct environmental threats, and thus it is difficult to predict the resistance functions of R genes based on synteny relationships to genes with known defense traits in other species. Discussion R genes from TD arrays are under more rapid evolution than R gene singletons. Whole genome synteny comparisons of R gene loci found that most R genes have syntenic orthologous R loci in other species. This suggests that R gene loci are relatively conserved among species, and R genes are mostly increased through TDs rather than by transposition or segmental duplications. These R gene TD arrays then evolved new functions against disease invasions. R genes featured in novel functions continuously evolve to adapt to changing environments and biotic disease challenges, so it is reasonable that R genes have evolved faster than other genes. However, when comparing the two groups of R genes (R TDs and singletons), we further found that syntenic families of R genes from TD arrays had higher K a values among genomes than R gene singletons. This suggests that R genes from TD arrays evolved faster than R singletons, and that TD arrays act as an important incubators for R genes to evolve temporary and variable functions to cope with novel disease factors. This conclusion is supported by the fact that many newly-mapped R genes are distributed in TD arrays, for example the blast resistance R gene Pi37 of O. sativa and the black streak disease resistance R gene Rox1 in Z. mays 24 . It is also consistent with previous findings that rapidly evolving R gene families contribute to the functional innovation of blast disease resistance after transforming them into rice 29 . K a /K s values further supported the conclusion that compared to R TD arrays, R singletons are under selection to be conserved. This indicates that when an important resistant-like function (long lasting, basic, or broad spectrum effect) has evolved, selection will retain the new variant in singleton status to prevent it from further sequence variation. These results support the hypothesis that the R genes from TD arrays and singletons adopt different strategies and play different roles in R gene evolution. R genes from TD arrays are apparently under diversifying selection in S. italica. The variome dataset generated from the resequencing of 916 accessions of S. italica 35 serves as a valuable resource for examining the features of genetic polymorphism in the S. italica population. By comparing the SNP and InDel variants between R genes and others, we found that R genes accumulated more mutations than other genes, which suggests that R genes evolved more rapidly than other genes. Furthermore, when taking these functional mutations into account, we found that R genes had significantly more such variants than non-R genes. Comparisons of R genes from TD arrays and singletons showed that R TDs had the highest average number of functional mutations per gene, suggesting that R genes from TD arrays are among the most changeable gene groups. Moreover, the differences in selection pressure imposed on R TD arrays, R singletons, and non-R genes resulted in R gene TDs having the highest level of polymorphism. For R singletons, there is selection for low levels of polymorphism, i.e., purifying selection. Considering that low-depth resequencing data cannot determine structural mutations, R genes may involve gene copy variants, despite SNP and InDel mutations. The current data may not reveal the full variation information of R genes in S. italica, but the data do reflect the general feature patterns. R genes are one of the most important gene families in plants, as they evolve relatively rapidly and adopt different strategies during the process of accumulating mutations and promoting functional innovations. In this work, we systematically examined the R genes in grass genomes and compared them across different species as well as in the population of S. italica using comparative genomics and population genomic analysis, respectively. We found that R gene loci show conservation among grass species, and that tandem duplication is the major pathway for R gene expansion in diploid genomes. More importantly, we found that R TD arrays and R singletons differed in their evolutionary strategy, with R TD arrays promoting functional innovations to cope with new disease threats, while R singletons were more likely to maintain and conserve important resistance functions over evolutionary time. Materials and Methods Resources of genome and gene sequences The 13 genome datasets (genome sequences, gene prediction coordinates, gene coding sequences and peptides) were downloaded from the databases phytozome ( https://phytozome.jgi.doe.gov/pz/portal.html ) and PlantGDB ( http://www.plantgdb.org/ ). The 11 grass species were: Triticum aestivum, Brachypodium distachyon, Brachypodium stacei, Oryza sativa, Oropetium thomaeum, Zea mays, Sorghum bicolor, Panicum virgatum, Panicum hallii, Setaria italica, Setaria viridis, and two species as outgroups, Musa acuminate and Arabidopsis thaliana. Resequencing data of 916 S. italica accessions were downloaded from the European Nucleotide Archive ( http://www.ebi.ac.uk/ena/ ) with the accession number ERP002070 35 . Determination of R genes Four main steps were applied to identify R genes in the 13 species. First, we retrieved the gene and protein sequences of 202 R genes from A. thaliana 8 and 456 R genes from O. sativa 22 and performed a Pfam_Scan to search for domains in these genes. Fifty-five domains were obtained. Second, we used a hidden Markov Model (HMM) search tool 34 to identify genes containing these 55 domains in the 13 species. Domains of identified genes were then extracted to check the domain compositions and orders in each gene. A subset of genes that contained the NBS domain were then chosen for further analysis. The N-terminal TIR domain was identified with the Pfam protein motif (PF01582; http://pfam.sanger.ac.uk/ ), while the N-terminal coiled-coil motif was identified by the tool NCOILS 39 . Tandem gene determination We searched all tandem gene arrays across each genome for sequence homology and locations of genes in local regions. The determining rule for a TD array was set as: 1) each tandem gene array is composed of continuously distributed homologous genes with a sequence identity of: Blastp E-value &lt; 1E-2), and 2) a tandem gene array should not be separated by insertion of more than five non-homologous genes. Synteny analysis of R genes Pairwise genome analyses were performed to detect orthologous gene pairs. First, we used Blastp to identify homologous genes (E-values &lt; 1E-20) between pairwise genomes. These homologous gene pairs were then submitted to check the homologous status of their flanking genes. Genes located in both flanking regions of the homologous gene pair were counted if they were the best hits of Blastp between the two species. Homologous gene pairs were considered as syntenic genes between the pairwise genomes if they had the highest ratio of homologous flanking genes. Alignment and phylogenetic tree analysis Protein sequences were aligned using MUSCLE 40 with default parameters. Then the aligned sequences were used to build a phylogenetic tree by MEGA using the Neighbor-Joining or Maximum-likelihood algorithms. Bootstrap tests were used to assign the confidence level of each branch of the tree. K a /K s analysis Protein sequences between pairs of homologous genes were aligned by MUSCLE 40 . Protein alignments were then translated into coding sequence alignments using an in-house Perl script. After that, K a and K s values were calculated based on the coding sequence alignments using the method of Nei and Gojobori as implemented in KaKs_calculator 41 . K a /K s values of syntenic orthologs between pairs of the 13 species were then plotted as histograms. Meanwhile, the K s loci were chosen from the aligned coding sequences and concatenated to build the phylogenetic tree. Variants calling and annotation S. italica resequencing data were analyzed following the method described previously 35 . Low-quality reads were removed: (1) if one read contained &gt;5% “N” bases; (2) if the average Phred-like score was &lt;20; (3) if less than 40-bp remained after a test trim of the 3ʹ nucleotides with a Phred-like score &lt;13; (4) if read-pairs were exact duplicates. Filtered reads were mapped to the reference genome of S. italica (version 1.0) using the “mem” algorithm of the Burrows-Wheeler Aligner (BWA) 42 . Variants were identified using Samtools 43 . Only loci with two alleles and Minor Allele Frequency (MAF) &lt; = 0.05 were kept for further analysis. For the information concerning gene prediction coordinates, variants were classified as being inter-genic, CDS, Intronic, or UTR. Variants in CDS were further separated into two types: those causing changes to coding proteins (non-synonymous SNPs and frame shift InDels) and synonymous SNPs and InDels without frame-shifts. Intronic variants were divided into splice site mutations (within 2-bp of splice site) and others. Genome diversity and selection We used two measures to estimate the selection pressure on the genome of S. italica populations: π and Tajima’s D. The statistic π measures genomic diversity 37 , computing the average difference per locus over each pair of accessions; Tajima’s D evaluates selection pressure to detect the genomic regions under purifying or balancing selection 38 . Both measures were calculated for variants in each S. italica gene.
Original work titled: 'Comparative Genomics Analysis in Grass Species Reveals Two Distinct Evolutionary Strategies Adopted by R Genes' (Scientific Reports/Nature Publishing Group UK) by the following authors: Yinan Zhang; Meijun Guo; Jie Shen; Xie Song; Shuqi Dong; Yinyuan Wen; Xiangyang Yuan; Pingyi Guo DOI: 10.1038/s41598-019-47121-8