The PanOryza pangene catalog of Asian cultivated rice

RPRP gene sets

We first collated gene models for the RPRP genomes as a source for the creation of a pan-Oryza gene set. All 16 genomes, including the IRGSP RefSeq were annotated using a consistent pipeline and data source (see Methods). In addition, there are also RAP-DB and MSU gene sets for the IRGSP RefSeq, making 18 annotation sets in total from which we built the pangene set. The counts of gene models for each genome used as input to the pangene set are provided in Table 1. In each accession, around 36,000 genes were identified except for Nipponbare, which has a greatly inflated gene count owing to the merging of three input sets (see “nipponbare_merged_GFF” at https://zenodo.org/records/14772953).

Pangene sets, clusters, and identifiers

Our first objective was to analyze gene annotation sets from the RPRP genomes and assign pangenes following WGA, indicative of synteny. The GET_PANGENES pipeline was used to create pangenes, determine their occupancy statistics, and assign long-term stable identifiers (Fig. 1). Supplemental Table S1 contains the matrix of pangenes, with their stable identifier (column 1) and then 16 further columns, one per input genome, containing transcript identifiers (if any) from each genome that have been mapped to that pangene. There are a total of 77,530 pangenes, including 35,029 singletons, that is, pangenes containing transcripts from only one genome (17,383 singletons when MSU gene models are removed, which heavily inflate the count), leaving 42,501 pangene clusters containing two or more members. Figure 1A displays the distribution of pangenes by occupancy class, demonstrating that most genes are core or cloud and that shell genes are relatively rare in the pangene set.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 1.

Pangene sets of rice. (A) Count of input genome members (x-axis) “occupancy” 1–16 versus counts of pangenes (pangenes containing only one MSU gene model have been removed). (B) Core-genome growth simulations after adding RPRP cultivars in random order (without singleton genes) using Tettelin (red) (Tettelin et al. 2005) and Willenbrock (blue) (Willenbrock et al. 2007) functions, fitted after 20 permutation experiments. (C) Counts of singletons, that is, genes per input genome present in clusters with no other gene models (for IRGSP, singleton count includes those with an OsNip [Gramene-generated] gene, for unbiased comparison). (D) Matrix of the percentage of shared clusters and heat map of RPRP pangenes. The dendrograms were computed by complete linkage clustering and Euclidean distances computed among columns.

We also assessed the quality and robustness of pangene formation through multiple metrics, including descriptive statistics on gene length, exon count, sequence distance between members, and “MSA completeness” (Supplemental Fig. S1; Supplemental Table S2), demonstrating that the majority of pangenes contain highly similar gene/protein sequences from different genomes. The metrics also allow for identification of pangenes containing one or more inconsistent members, indicative of different overlapping gene predictions in some genomes (e.g., coding sequences predicted on different frames), which can be artifacts of incorrect gene model prediction.

Figure 1B shows that the core set of genes converge at 25,178 genes according to the Tettelin method (23,987 by Willenbrock method), indicating this is the minimal core set contained within an O. sativa genome. The pangene curve shows that the pangenome size starts converging at around 42,000 genes (Supplemental Fig. S2). Figure 1C displays the counts per variety of singletons from each genome, ranging from 457 to 1511 (mean, 766; median, 752), indicating consistency across different cultivars for the presence of these unique genes. For Nipponbare, the results are filtered to include singletons containing an OsNip-source model, created using the same approach as the other 15 genomes to avoid artefactually high counts owing to merging different input gene sets, giving a set of 964 genes in line with other genomes. Figure 1D shows the pangene cluster overlap across cultivars. In general, there is reasonable clustering by rice subfamily (japonica vs. aus vs. indica). The merged Nipponbare annotation has lower overlap with other clusters, owing to it having a much higher count of input genes.

Genomic positions of pangenes

The GET_PANGENES approach uses WGA and does not enforce that genes must be located on the same chromosome. This enables the algorithm to identify cases in which breeding (or natural selection) has caused regions of the chromosome to translocate or potentially in which there have been assembly errors. In Supplemental Figure S3, we display the chromosomal location of genes within the pangene set against the Nipponbare (IRGSP) reference genome (see Methods). There are sporadic cases of genes that are not collinear with respect to the reference genome, with a few notable systematic cases. These include a set of genes that are inverted on Chromosome 1 for OsLima at position 600,000–800,000 bp with respect to all other genomes. Similarly, there are two inversions on aus genomes (N22 and Natal Boro) around the center of Chromosome 1. Another notable case is present on Chromosome 6, in which Nipponbare has a large inversion compared with the other 15 genomes. Data containing the positions of genes within pangene clusters can be found in Supplemental Table S3. There are 1502 pangenes that have occupancy greater than two and contain genes mapped to different chromosomes (Supplemental Table S3), which researchers should consider when using genomics resources, for example, SNP mapping or allele mining.

Pangene distribution across rice families

Supplemental Table S4 contains the proportions of genomes by subfamily (japonica, aromatic, aus, indica) that have a member within each pangene. For example, by filtering these data, one can identify genes that segregate by family, for example, 359 pangenes found in all japonica genomes but only zero or one other genome, 229 pangenes found in eight or nine indica genomes but present in zero or one other genomes, and 336 pangenes found in aus varieties and zero or one other genomes. For multiple research and breeding purposes, such genes are of high interest, particularly if there are QTLs mapped to these genomic regions. One of the current major limitations of genome-wide association or other population SNP mapping studies is that the ISRFG RefSeq reference genome (IRGSP) is almost always used for read mapping and locating genes near SNPs. Unfortunately, as we show here, there are 17,130 pangenes that do not have an IRGSP gene model (5634 with occupancy of two or more), indicative of a large pool of genetic resources that is currently missed in most such studies.

Properties of pangenes by occupancy

We next explored the properties of pangenes related to occupancy (Fig. 2). We chose to focus on quality metrics derived from the RAP-DB reference gene annotations (for which most data already exist from different resources) and take the occupancy statistic from the parent pangene. First, we assessed the quality of a given gene model, using the TRaCE (Olson and Ware 2021) algorithm, and assigned a minimum annotation edit distance (AED) score, from three transcriptomic data sets, to assess the quality of the annotation. AED is calculated on a scale of zero to one based on the sequence distance between the gene model and the best transcript, assembled from publicly available transcriptomics data.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 2.

Exploration of pangenes and their properties by their occupancy. The x-axis shows the pangene occupancy for each cluster; the y-axis shows data derived from RAP-DB gene models within those clusters (clusters containing no RAP-DB model are excluded). (A) Mean of minimum AED score from three transcriptomic data sets. (B) Proportion of genes that contain a significant match to an InterPro domain. (C) Mean count of amino acids. (D) Boxplot showing the counts of orthologs derived from Ensembl Compara (see Methods). (E) Boxplot showing counts of paralogs. (F) Boxplot of mean quantitative gene expression data log₂(FPKM + 1), sourced from 11,726 rice RNA-seq samples. (G) Boxplot of mean prediction score (pLDDT) from the corresponding AlphaFold2 model. (H) Percentage of clusters with genes having peptide evidence (n > 2 per gene). The p-value indicates result from Tukey's HSD test for core versus cloud pangenes.

Figure 2A shows the mean (minimum) AED score for all RAP-DB genes, against the occupancy of their parent pangene. Core (occupancy, 16) or softcore (15 or 16) clusters have significantly lower AED within core (mean, 0.5) versus cloud pangenes (mean, 0.9), indicating that core genes are called with higher confidence using transcriptomic data sources. Of the core pangene, 83% have an assigned InterPro domain (Fig. 2B), compared with 11% in singleton pangenes (occupancy, 1). A more detailed analysis of protein domains across all pangenes is provided below. Core pangenes are also significantly longer and are more likely conserved across the plant kingdom as indicated by their protein length and ortholog count (Fig. 2C,D).

There is a more complex trend for paralog count (Fig. 2E), in which there is no significant difference between paralog counts for cloud or shell versus core sets. There appears to be a trend in which shell genes have higher paralog counts. We hypothesize that there are two “competing” distributions here. We would expect that pangenes of lower occupancy are more likely to have paralogs, as this set contains rapidly duplicating gene families (further discussed below), but are also enriched for faulty gene models in which paralogs cannot be identified.

Pangenes with softcore/core occupancy have significantly higher gene expression values compared with cloud genes, with singleton pangenes (occupancy, 1) frequently having no detectable signals; these could be pseudogenes or genes only transcribed under particular conditions not assessed in the source data (Fig. 2F). Prediction of a 3D protein structure for core/softcore and shell pangenes can be done by AlphaFold with high confidence (median ± interquartile range: 78.1 ± 20.6 for core, 75.0 ± 23.4 for softcore, and 70.3 ± 27.2 for shell pangenes), whereas a protein structure for cloud pangenes has a median score of 54.7 ± 16.6, indicating lower confidence (Fig. 2G). Cloud genes thus appear to be enriched with low-quality gene models (short sequences, no InterPro domains, weaker structure models), which is further supported by peptide-based evidence derived from experimental proteomics studies (further details below), in which ∼9.9% of singleton pangenes are supported by more than two peptides per gene. In contrast, 52.1% core RAP-DB genes show more than two peptides per gene as evidenced for the pangene cluster (Fig. 2H).

Support for rice panproteome

To determine experimental support for the rice pangenes, we performed large-scale reanalysis of public proteomics (mass spectrometry) data, against a comprehensive protein database derived from all possible gene models within the pangene set. We searched 19 public data sets (about 30 million MS/MS spectra) sourced and pooled from multiple rice varieties to maximize total coverage. All data sets were from “shotgun” methods; that is, proteins were digested into a total peptide pool prior to LC-MS/MS. As a result, for peptides that can match to more than one gene product (per genome), it is not straightforward to determine which is the correct assignment. Using an all peptide-to-protein mapping, we demonstrate potential peptide-based evidence using more than 200,000 distinct peptides for 435,196 gene models and 293,844 genes present in 28,658 pangenes.

Taking a parsimonious approach in which a peptide can only be mapped to one protein, we found around 13,000 genes in each genome with at least two peptides, indicating strong evidence for the gene's protein-coding potential (Fig. 3A). In agreement with the previous set of observations, we find a significantly higher number of peptides mapped to core and soft-core pangenes than to cloud pangenes (Fig. 3B). Out of all the peptides mapped to all the proteins, ∼73% (n = 168,151) of the peptides mapped to proteins in all 16 rice accessions, indicating good support for these pangenes, and in general, based on current evidence, there is little deviation in sequences for readily detectable rice proteins across the pangenome (Fig. 3C).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 3.

Experimental proteomics data support for pangenes. (A) Count of proteins per annotation set with at least two supporting peptides; (B) log counts of peptide sequences versus genome occupancy; (C) absolute counts of peptides mapped to proteins from one to 16 genomes; and (D) scatterplot of mean (min) AED score across three transcriptomics data sets for those (RAP-DB) gene products with or without peptide evidence, by genome occupancy of the pangenes. Dashed line indicates linear regression, indicating that gene products with peptide evidence have better (lower AED) support from transcriptomics data. The p-value indicates result from Tukey's HSD test for core versus cloud pangenes.

We further demonstrate the relationship between proteomic evidence (those genes with or without proteomics evidence) and AED score (derived from transcriptomes mapped to the RAP-DB gene representative of each pangene cluster) in Figure 3D. Proteins without experimental support tend to have high AED scores (ranging from one to 0.8, depending on genomes occupancy), whereas those with experimental support and high occupancy have low AED scores, indicative of well-supported, high-confidence gene models.

Exploring the predicted protein domains of pangenes

We ran InterProScan (Jones et al. 2014) over all input proteins available within the pangene set, assigning Pfam (Mistry et al. 2021) and InterPro (Paysan-Lafosse et al. 2023) domains. For each pangene, we assigned additional data types and statistics, including the total number of different domains identified, the genomes (and proteins) containing those domains, and a representative domain for the cluster (most common across genomes), available in Supplemental Table S5. We then assigned which proteins within a pangene cluster contain the representative domain, and which proteins are outliers, that is, do not contain the most common domain. Of the 42,501 pangenes with at least two members, 26,044 have at least one assigned Pfam domain (and 28,507 at least one InterPro domain). Of these, 21,779 are “domain-consistent” clusters (83%), having zero genomes lacking the most common Pfam domain. Similarly, 88% of pangene clusters have a consistent common InterPro domain. In Supplemental Table S6, we show the counts of proteins that do not contain the most common InterPro domain per cluster per input genome for the 3485 pangenes that are not “domain consistent.” There are around 400–600 proteins per input annotation set that do not match the most common domain. RAP-DB and MSU gene models have higher counts but have not been annotated using the same pipeline as other genomes and thus are less consistent than OsNip models compared with other genomes, and we cannot straightforwardly conclude that the gene models from these sources are less reliable.

We identified the most enriched domains found in the set of core versus cloud pangenes and observed some distinct patterns of domain types (Fig. 4A). Core pangenes are enriched in domains including S-adenosyl-L-methionine-dependent (SAM) methyltransferase, AP2/ERF, WD40, helix–loop–helix, Myb domain, C2H2 zinc finger, and homeobox domain. Although SAM methyltransferases are responsible for methylation of several substrates, other domains among these pangenes are commonly present in DNA-binding transcription factor (TF) proteins. This enrichment among core pangenes thus reflects the conserved nature of these molecular processes. However, the domains within the cloud pangenes appear to be more varied with transposases seen as a prominent group. The enrichment among cloud pangenes may be indicative of genes that have either a nonconserved nature, or they failed to cluster for a variety of technical reasons.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 4.

Exploration of predicted pangene domains. (A) Enriched InterPro terms among cloud versus core pangenes. The top 30 terms are shown in the plot. Size of the dot represents the enrichment factor, and the color of the dot shows the P-value. (B) Histogram shows the counts of domains within the pangenes versus their mean domain occupancy for both the Pfam and InterPro domains. The domains were filtered for presence in a minimum of five pangenes. (C) The plot shows the percentage of pangene clusters according to their occupancy status for two gene families: NAC domain-containing TFs and NB-ARC domain-containing genes. (D) The variation in mean domain occupancy of various Pfam domains for the selected major gene groups is shown. The solid back dot shows the median value. (E) The scatterplot shows the mean domain occupancy of various domains under the group “domains with unknown function” (DUFs) versus the number of pangenes with the respective domain. DUFs present in more than 50 pangenes are highlighted in black text.

Among the significantly enriched domains in the cloud pangenes, we found a domain named DNA-directed RNA polymerase III subunit RPC4 (PF05132/IPR007811) that is involved in the transcription of small RNAs (Fig. 4A). There are only four RPC4 pangenes with occupancy = 16, a varying number of shell pangenes, and 228 cloud RPC4 pangenes with occupancy = 1. In Supplemental Figure S4, we show the expansion of this gene family in indica- and aus-type rice varieties with 15–38 pangenes, whereas japonica/aromatic rice varieties have only two to three pangenes of this domain, potentially indicating a very substantial expansion of this family in indica- and aus-type varieties.

Next, we explored the occupancy of protein domains across the pangene set. Variance in the distribution of domains can potentially inform on functional units that either have a tendency to become duplicated or have been more frequently introduced via introgression, compared with those that are invariable. Figure 4B contains a histogram of the mean occupancy of protein domains (InterPro and Pfam) showing a multimodal distribution. As noted in the Methods, we then split the distribution into three groups: “invariable (occupancy) domains,” “partially variable domains,” and “highly variable domains.” We hypothesize that gaining (or losing) copies of invariable domains is highly detrimental to rice. GO enrichment analysis of pangenes containing these groupings of domains provides support for this hypothesis (Supplemental Fig. S5). For example, pangenes containing ubiquitin-protein transferase and TF DNA binding activity show conserved and invariable domains, whereas domains involved in cell surface receptor signaling appear highly variable, indicating the fast-evolving nature of these domains driven by selective breeding pressure, like disease resistance.

As a case study to test the pangene clusters, we analyzed two gene families (Fig. 4C): the NAC TF family (an invariable domain) and the nucleotide-binding domain shared by Apaf-1, certain R-proteins, and CED-4 (NB-ARC) disease-resistance family (a highly variable domain) (Jacquemin et al. 2014). NAC TF family genes have been characterized for their roles in abiotic stress responses in rice (Nakashima et al. 2012). The NB-ARC domain is a conserved nucleotide-binding domain found in a variety of proteins across different organisms, including plants, animals, and fungi. It plays a crucial role in processes like apoptosis (programmed cell death) and innate immune responses. We identified 175 NAC pangenes with 142–154 genes in each genome, and the domain has a mean occupancy of 13.2. 77% of NAC pangenes (n = 135) have core occupancy and 9% cloud (n = 16), indicating a highly invariable domain (Fig. 4C; Supplemental Fig. S6). In comparison, the NB-ARC is a highly variable domain (mean occupancy, 7.2) identifiable within 1005 pangenes with 432–518 genes per genome. NB-ARC family genes are mostly found with an occupancy of 16 or one, indicative of rapid introduction of additional copies or recent gene duplications in some genomes (Fig. 4C; Supplemental Fig. S7). There is also a difference between NAC and NB-ARC in shell pangenes: 9% of NAC pangenes (n = 16) versus 34% of NB-ARC genes (n = 341). One could hypothesize that because NB-ARCs have a role in pathogen effector recognition, it is desirable for the plant to have a wide diversity of genes but to tolerate the loss of a small number of genes in some genomes (leading to shell pangenes), whereas loss of a TF is highly detrimental, and the existence of cloud genes is explained by rare duplication events.

A similar trend is seen with other domains like leucine-rich repeats (LRR; median occupancy, 7.50 ± 0.945), indicating the highly variable nature of these domains (Fig. 4D). Given the role of LRRs in plant immune responses, the domain occupancy data are suggestive that these genes are fast evolving through breeding, for example, when past gene duplications giving resistance to a given pathogen in cultivated rice have been selected. Notably, it is one of the domains associated with pangenes enriched with the highly variable GO term “cell surface receptor signaling pathway” (Supplemental Fig. S5). Multiple domain types with the name “F-box” show a partially variable occupancy (9.87 ± 2.66) compared with the least variable domains of TF (13.3 ± 1.73) or protein kinase (13.4 ± 3.70) types. F-box domain-containing proteins include those involved in the SCF type ubiquitin E3 ligase family and, in plants, have multiple functional roles, including stress responses. Our data also suggest past selective breeding pressure, giving rise to variability in the domain occupancy across the rice pangenome.

The class of “domains with unknown function” (DUFs) shows a median occupancy of 13.6 but with a wider interquartile range (±5.49), with 46 domains present in more than 10 pangenes (Fig. 4E). Several DUFs (DUF594, DUF4220, DUF295, DUF1618, DUF6598, and DUF1668) are present in more than 100 pangenes and have been classified as highly variable domains, which we believe may be indicative of past selection. The genes containing these domains are thus candidates for functional validation and are potential carriers of desirable traits.

Infrastructure for pangene exploration

The pangene set has been loaded into popular, widely used databases, including Ensembl Plants and Gramene, and into our own pangene resource (https://panoryza.org/). At Gramene (https://oryza.gramene.org/), it is possible to search with pangene identifiers, for example, Os4530.POR.1.pan0020022, which is the pangene identifier for NAC106 (RAP-DB: Os08g0433500) (Fig. 5A). NAC106 is a NAC TF, which has been associated with salt tolerance, tiller angle, and leaf senescence (Sakuraba et al. 2015). The Gramene search then returns the syntenic and orthologous gene records within each of the genomes assigned to that pangene identifier. There are multiple views for exploring the data, including as phylogenetic trees within the set or as expansions include orthologs from other plants. The RPRP assemblies and annotations have been added into Ensembl Plants, and the new Ensembl site which is currently in beta (beta.ensembl.org). The pangene identifiers are visible as gene synonyms on beta and will become available via Ensembl Plants in a future release (Fig. 5B).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 5.

Representative example of infrastructure for pangene exploration using pangene Os4530.POR.1.pan0020022. (A) The screen shot shows the record of pangene Os4530.POR.1.pan0020022 on the Gramene database (https://oryza.gramene.org/). (B) The page shows the gene record for a member gene OsARC_08g0017270 of Os4530.POR.1.pan0020022 on the Ensembl beta release.

At the PanOryza site, we have a JBrowse (Buels et al. 2016) instance in which each genome can be explored with various tracks of data (Supplemental Fig. S8A). For the IRGSP reference, we have aligned gene models from other cultivars using LiftOff, enabling any variation in models from other genomes to be displayed. For NAC106 as an example, which has occupancy = 16, some variation is suggested within exon 2 in the N22, IR64, and Azucena varieties. Another example is shown in Supplemental Figure S9 that shows the Gramene search result and PanOryza genome browser (JBrowse) for Os4530.POR.1.pan0035634. The members of this pangene coding for NB-ARC domain-containing proteins are present in all japonica and Basmati accessions of RPRP but absent in aus or indica genomes.

We have also created an R library, with code for interactive exploration of clusters using R Shiny through a web browser. The application allows users to visualize the relative position of genes within a cluster across the 16 genomes and explore the members and metadata about each pangene (Supplemental Fig. S8B; https://github.com/PGB-LIV/PanOryza-pan-genes-release-v1.0/tree/main/heatmap_app/).

A permanent DOI has also been created for the pangene matrix and identifiers (https://zenodo.org/records/14772953) so that they can be used now in other research projects. It is anticipated that they will be updated in the future (e.g., yearly), under some of the following circumstances: (1) gene models change; (2) additional genomes are added; and (3) there are refinements to the algorithm for creating pangenes. Full guidelines for the mechanism of updating pangene identifiers are also available at the Ensembl Plants GitHub repository (https://github.com/Ensembl/plant-scripts/tree/master/pangenes). Mappings will be made available so that any users can track how and why pangenes have changed. These resources can all be used by the rice research and breeding community now to assist in mining the wide pool of genetic resources available for Asian rice.