Pan-genome inversion index reveals evolutionary insights into the subpopulation structure of Asian rice

The 18-genome data package

To investigate the genome inversion landscape of Asian rice from a population structure perspective, we first combined a set of 16 previously published high-quality genomes32,33,34 that represent the K = 15 population structure of O. sativa, plus the largest Xian/indica (XI) admixed subpopulation (XI-adm: Minghui 63 (MH63)) to create a “Rice Population Reference Panel” (RPRP – yongzhou2019.github.io/Rice-Population-Reference-Panel/). This panel was annotated with a uniform pipeline to minimize methodological artifacts using RNA-Seq and Iso-Seq data generated or collected from all 16 genomes as evidence (Table 1, Supplementary Tables 13, Supplementary Fig. 1, Supplementary Data 1, and Supplementary Note 1).

Table 1 Assembly and annotation statistics of the 18-genome data package

To anchor this panel within a phylogenetic context, we long-read sequenced, and de novo assembled two additional genomes from both a representative species of the progenitor of Asian rice – i.e., O. rufipogon (AA) and the African BB genome outgroup species – O. punctata (Table 1, Supplementary Table 4, and Supplementary Fig. 2). Both species are diploid and have similar genome sizes as Asian rice. These genomes were assembled to a similar quality as the RPRP (i.e., BUSCO > 95%, number of gaps <50, ContigN50 > 10 Mb) (Supplementary Table 4 and Supplementary Note 1).

All 18 genomes and their annotations are henceforth referred to as the “18-genome data package” (See “18-genome data package” in the Supplementary Note 1 section for a complete description of this data set).

Creation of a pan-genome inversion index for Asian rice

To ensure the detection of the majority of inversions present in cultivated Asian rice, we selected an additional 57 publicly available long-read genome assemblies (out of 94) for analyses that met similar quality standards present in the 18-genome data package, for a total of 75 genomes (Supplementary Fig. 3, Supplementary Data 2, and Supplementary Note 2).

To detect inversions, we compared pairwise all 74 reference genome assemblies with the IRGSP RefSeq34, following a workflow as detailed in Supplementary Note 3 (Supplementary Fig. 4 and Supplementary Table 5), and identified a total of 12,141 inversions (≥100 bp) (Supplementary Data 3), 1769 of which were non-redundant (Supplementary Data 4). Finally, since several inversions had overlapping coordinates (i.e., 80% of their length), they were clustered to obtain a set of 1426 clustered inversions (Supplementary Data 4). We then validated a subset of 264 inversions from four randomly selected genome assemblies (i.e., LM, NABO, CM, and MH63) using PacBio long-reads, and found that 97.7% of the inversion breakpoints could be detected (Supplementary Table 6), thereby confirming the high accuracy of our detection strategy.

To determine if the use of 75 genomes was sufficient to call the majority of 100 bp or greater inversions in the pan-genome of Asian rice, we performed a permutation test (n = 1000), and found that the use of ~60 genomes was sufficient to capture the majority of inversions with allele frequencies greater than 2 out of 75 genomes (Fig. 1a, b). These results demonstrate that employing 75 genomes, that bridge the K = 15 population structure of cultivated rice, can yield a robust pan-genome inversion index for Asian rice.

Fig. 1: Rice inversion index summary.
figure 1

a, b Resampling permutation test to identify the relationship between the number of genomes and all inversions and shared (non-genome specific) inversions, respectively; c Density of inversion lengths; d Bionano validation of inversions larger than 1 Mb, i.e., Clu-INV0100180, Clu-INV0100660, and Clu-INV0600550. In each panel, the top line shows the optical map used as a reference, the bottom line shows the genome assembly of the variety with the inversion. Gray lines connect restriction sites that are aligned (blue regions), while yellow segments show unaligned regions. Black boxes highlight the position each inversion. Source data are provided as a Source Data file.

Our inversion index was compared with previous studies in rice, i.e., 3K-RGP11 and the pan-genome study of 33 rice genomes12. Inversions were treated as identical if they matched the following two criteria: 1) the inversion length difference was smaller than 200 bp, and 2) the differences between the coordinates of breakpoints of inversions were smaller than 100 bp. In doing so, we found that of the 1769 nonredundant inversions identified, 38.6% have already been reported (Supplementary Data 4).

As expected, more inversions were observed when we compared the BB genome outgroup (O. punctata) than AA genomes to the IRGSP RefSeq: 337 (total length = 15.16 Mb) (Supplementary Data 3), followed by the AA genome outgroup (O. rufipogon): 230 (total length = 12.95 Mb). On average, each O. sativa genome was found to contain 164 inversions, ranging from 33 (TG22) to 230 (YX1) (Supplementary Data 3). We found a larger number of inversions (i.e.,164–230) when comparing the O. sativa XI-subgroup genomes to the GJ-temp: IRGSP RefSeq than when comparing the O. sativa GJ-subgroup genomes (33–201 inversions) to the same reference, which is consistent with the long divergence time estimate between the XI and GJ varietal groups (i.e., 360,000–400,000 years35).

The total length of inversions from one genome to another ranged from 14.95 Mb (cA2: NATEL BORO) to 255 Kb (GJ-temp: TG19) (Supplementary Data 3) and totaled 117.75 Mb across all 75 genomes tested. We identified 470 clustered inversions shorter than 1 Kb and 11 clustered inversions greater than 1 Mb relative to the IRGSP RefSeq (Fig. 1c and Supplementary Data 4). These extra-large inversions (>1 Mb) were found on seven chromosomes (i.e., 1, 4, 5, 6, 8, 10, and 11) and included seven of which were previously reported (i.e., Clu-INV0600550 (4.5 Mb)30, Clu-INV0800800 (1.1 Mb)36, Clu-INV0100660 (1.8 Mb)12, Clu-INV0500270 (3.1 Mb)12, Clu-INV0800450 (5.2 Mb)12, Clu-INV0800470 (1.2 Mb)12, Clu-INV110070 (1.1 Mb)12, and four were identified in this study (i.e., Clu-INV0100180 (2 Mb), Clu-INV0400050 (1.3 Mb), Clu-INV0400650 (4.3 Mb), and Clu-INV1000940 (1.3 Mb)). Three of the 11 clustered inversions (i.e., Clu-INV0100180, Clu-INV0100660, and INV0600550) were validated with available Bionano optical maps (Fig. 1d), while the remaining eight were not due to the lack of Bionano data.

Chromosomal distribution of the pan-genome inversions index

To determine if the inversions detected were uniformly distributed across all chromosomes, we tested for uniformity using the Kolmogorov–Smirnov (KS) test, and found that only one chromosome, chromosome 4, deviated significantly from uniformity (i.e., p = 3e−04; Fig. 2a, b, and Supplementary Data 5). To check if inversion size played a role in this chromosome 4 anomaly, we repeated the KS test by classifying inversions according to their lengths (<1 Kb, 1–5 Kb, 5–10 Kb, and >10 Kb). We found no significant deviation from uniformity except for the >10 Kb inversion class on chromosome 4 (Supplementary Fig. 5 and Supplementary Data 5), where inversions are mainly concentrated near the centromere on both the short and long arms of chromosome 4. These results demonstrate that the inversions detected appear evenly distributed genome-wide, with one exception.

Fig. 2: Genome-wide inversion distribution.
figure 2

a Chromosome distribution of the pan-genome inversion index, and inversion hotspots. Chromosome heatmaps in the middle along the chromosome represent the density of inversions. The line on the left of each chromosome represents the number of inversions per 200 Kb window. The red dotted line cut off is the top 2% number of inversions in per 200 Kb window. The blue boxes on the right of each chromosome represents the inversion hotspot regions. b The Kolmogorov–Smirnov (KS) test for inversion uniformity distribution across the 12 rice chromosomes. The black line is actual inversion distribution, and the red line is the 10,000th uniformly distributed simulation. Source data are provided as a Source Data file.

To search for inversion hotspots, we performed a 200 Kb window-based analysis across all chromosomes. We defined the top 2% of all windows with the highest frequency of inversion start coordinates as hotspots. This analysis revealed 47 putative hotspots, including 239 independent inversions where 12 inversions overlapped with centromeres (on chromosomes 7 and 8) but none with telomeres (Fig. 2a and Supplementary Data 6).

Phylogenetic analysis of the pan-genome inversion index

To determine if the pan-genome inversion index could be used to phylogenetically distinguish each of the 75 high-quality genomes into the expected K = 15 subpopulation structure, we used the unweighted pair group with arithmetic mean distance tree method, with the O. punctata (BB) and O. rufipogon (AA) genomes as outgroups (Fig. 3). As shown in Fig. 3, 66 of the 73 genomes could be subdivided across the expected subpopulation structure of Asian rice, while the remaining seven fell into two XI clusters that have yet to be characterized.

Fig. 3: Phylogenetic tree of 75 high-quality rice genomes using the pan-genome inversion index.
figure 3

Phylogenetic relationships of the 75 high-quality genomes used to create the pan-genome inversion index, inferred using the UPGMA method (unweighted pair group method with arithmetic mean95). The tree is drawn to scale, with branch lengths in the same units as those of the evolutionary distances used to infer the phylogenetic tree. Evolutionary analyses were conducted in MEGA796. Two wild relatives (O. rufipogon and O. punctata) were heighted with light gray arc. Subpopulations from GJ, cB, XI, and cA groups were highlighted with light blue, dark blue, light yellow and dark yellow, respectively. Two uncharacterized XI subpopulations are shown with a black circle. The distance matrix of SNP and INV polymorphisms was significantly correlated (Mantel test, simulation with n = 999,999, r = 0.79, p = 1e−6). Source data are provided as a Source Data file.

Inversion rate estimations for Asian rice

It was estimated that the AA genomes of the Oryza diverged from the BB genome type about 2.5 million years ago (MYA)2, which equates to an inversion rate of 67.4 inversions per million years (Table 2 and Fig. 4a). If we use the inferred AA genome diversification (i.e., O. sativa GJ-temp IRGSP vs. O. rufipogon) rate of ~0.50 net new species/million years2,31, then the estimate is 230 inversions per MY (Table 2 and Fig. 4a).

Table 2 Inversion rate estimations
Fig. 4: Species-specific, population-specific and shared inversion analysis of the pan-genome inversion index of rice.
figure 4

a A model showing species specific inversions and inversion rates in Asian rice and two wild relatives (O. punctata and O. rufipogon). b Frequency of genome-specific, subpopulation specific, group specific and group shared inversions across the population structure of Asian rice (n = 631). Source data are provided as a Source Data file.

Regarding O. sativa subpopulations, it was estimated that temperate japonica (GJ-temp) first diverged from proto-japonica about 14,200 years ago37, thus we used this population divergence time for comparisons within Geng/japonica. For the computation of time to a MRCA for a pair of genomes, we added the expected time to coalesce within the proto-japonica population38, which can be estimated to be 10,000 generations based on previous estimates of effective population size37. For a comparison involving tropical japonica, care should be taken to exclude inversions that possibly originated in Xian/Indica through introgression. Thus, out of the 431 clustered inversions that were segregated within GJ, we focused on 170 that were not segregated within non-GJ subpopulations. To deal with the possibility of inversions that originated in XI but did not appear in our data set, we also filtered out 26 inverted regions closer to XI than to GJ based on the 3K-RGP SNP data39. The remaining 144 inversions were used to compute the number of inversions between each pair of GJ genomes. The average number of (filtered) inversions between two GJ genomes is 20.86, which translates to an inversion rate of 735 inversions per MY. Calculations involving comparisons only between a temperate GJ and another GJ population lead to a similar estimate of 749 inversions per MY (Table 2 and Fig. 4a).

Species and population scale analysis

Of the 1426 clustered inversions detected (Supplementary Data 4), we classified them into four different groups: Inversions segregating in O. sativa (S), Inversions segregating in both O. sativa and O. rufipogon (SR), O. rufipogon specific (R), and O. punctata specific or AA-fixed (i.e., ancestral state not clear) (P) (Fig. 4a). As a result, 1303 (91.4%) appeared to be species-specific, i.e., O. sativa (S): 885 (totaling 54.07 Mb), O. rufipogon (R): 96 (totaling 4.66 Mb), and O. punctata specific or AA-fixed (P): 322 (totaling 10.06 Mb) (Fig. 4a). The O. sativa specific inversions (S) were further classified into three categories: i.e., the IRGSP-1.0 RefSeq represents ancestral state (S1), IRGSP-1.0 RefSeq represents the derived state (S2), all O. sativa represent the derived state (fixed inversions) (S3), in which 872 (48.7 Mb), 11 (5.2 Mb), and 2 (72.5 Kb) inversions for S1, S2 and S3 were observed, respectively (Supplementary Table 7). The remaining 123 non-specific inversions were found to segregate in both O. sativa and O. rufipogon (i.e., originated before O. sativa and O. rufipogon divergence, or introgression) and totaled to about 3.3 Mb in size, including 121 inversions with IRGSP-1.0 RefSeq as the ancestral state (SR1), and 2 inversions with IRGSP-1.0 RefSeq as the derived state (SR2) (Supplementary Table 7).

To gain deeper insight into the distribution of inversions in rice at the population scale, we studied the 872 O. sativa specific inversions (IRGSP-1.0 RefSeq has ancestral state) across the 3K-RGP data set28. Beginning with a manually curated set of 8720 inversion alignment patterns at their breakpoints (i.e. 872 inversion alignment patterns across ten accessions), we classified their short-read alignment patterns into four categories (i.e., presence of both breakpoints = inversion, absence of both breakpoints = no inversion, presence of a single breakpoint = inversion + deletion, and no data at both breakpoints = NA) (Supplementary Fig. 6a). These categories were then used to train a machine learning workflow to search for the presence or absence of 2,636,928 alignment patterns for 872 inversions across the 3K-RGP data set (Supplementary Fig. 6b, c and Supplementary Table 8). Of note, this analysis identified 241 inversions and 273 accessions with more than 30% missing inversion data and were thus filtered out, leaving a final data set of 631 inversions across 2751 accessions for downstream analysis (Supplementary Data 7 and Supplementary Note 4). These 631 inversions were then classified into genome-specific, subpopulation-specific, group-specific, and group-shared inversions based on their inversion frequencies (Supplementary Data 8). As shown in Fig. 4b, 94 genome-specific, 49 subpopulation-specific, and 37 group-specific inversions (including 30, 3, and 4 specific inversions to XI, cA, and GJ groups, respectively) were observed. We also observed that 451 inversions were shared among different groups, of which 386 (61.1% of total O. sativa specific inversions) were shared across the GJ, XI, cA, or cB groups, 58 were shared between the XI and cA groups, and seven were shared between the GJ and cB groups (Fig. 4b). In total, these results reveal that 91% of the 1426 cluster inversions tested were species-specific (i.e., O. punctata, O. rufipogon, and O. sativa) at the population scale, and may provide clues to their possible roles in speciation over evolutionary time.

Characterization of transposable element content within inversions and their breakpoints

Transposable elements (TEs) are known to be associated with inversions11,12. Thus, we analyzed the TE content across the inversion index and at their breakpoints. The total amount of TE related sequences present in the pan-genome inversion index was 63.4%, which is significantly higher (Student’s test, p = 0.0003) than the average TE content present in the 18-data genome data package at 51.3% (Supplementary Data 9). Furthermore, out of the complete set of 1769 inversions, 888 showed partial or total similarity to TEs at both breakpoint regions (+/−100 bp from the breakpoints), and another 389 showed partial or total similarity to TEs in at least one breakpoint. These results demonstrate that TEs are enriched within inversions and their breakpoints in Asian rice.

Analysis of these breakpoints revealed that both long terminal repeat retrotransposons (LTR-RTs, i.e., Ty3-Gypsy and Ty1-Copia) and DNA TE Mutator-like elements (MULEs) were significantly enriched (student’s test, p = 5.08E−11 to 0.027), when the frequency of their presence at the 3538 breakpoints was compared to 35,380 randomly selected genomic locations (i.e., ten replicates) (Fig. 5a). We further studied TEs at the breakpoint of each inversion that were shared across all Asian rice genomes. In doing so, we identified 17 TE families (i.e., 13 Ty3-Gypsy, 1 Ty1-Copia, 2 CACTA, and 1 Mutator) present at the breakpoints of more than 10 inversions (Fig. 5b and Supplementary Data 10). An example of an inversion enriched in TEs, including the internal and LTR portions of at least three different LTR-RTs, is shown in Fig. 5c.

Fig. 5: Transposable element analysis of the pan-genome inversions index of rice.
figure 5

a The amount (y-axis) of different transposable element (TE) families (x-axis) show that three TE families (i.e., LTR-RT Ty1-copia, Ty3-gypsy and DNA-TE MULE) were observed in higher frequencies at the breakpoints of the pan-genome inversion index than the resampled control tests. Box- and bar-plots show the frequencies of TEs observed at the breakpoints of random resampled regions (n = 10 biologically independent resamples) and the pan-genome inversion index, respectively. Each boxplot presents the minimum, first quartile, median, third quartile, and maximum value, and along with mean ± SD (Standard Deviation) are shown. b Enrichment/depletion of 17 TEs present at the inversion breakpoints with more than 10 copies. c Details of Ty3-gypsy Os0025 presence at inversion breakpoints, with support from PacBio long reads. Os0025_LTR (long-terminal repeat), and Os0025_INT (integrase). Source data are provided as a Source Data file.

Since it is known that inverted repeats can trigger ectopic recombination, thereby leading to genome inversions40,41, we interrogated 100 randomly selected O. sativa inversions for direct or inverted repeats in close proximity to inversion breakpoints. This analysis revealed the presence of direct repeats (including microsatellites) at the breakpoints of both ends of 11 inversions, and inverted repeats at single breakpoints of 30 inversions (Supplementary Data 11). For the remaining 118 inversion breakpoints (i.e., 59 inversions), no clear evidence could be found for the presence of inverted repeats at their inversion breakpoints. Of note, our analysis only considered inverted repeats longer than 10 bp, even though shorter have been shown to trigger inversions as well42. For this reason, we caution that our estimate of the presence of inverted repeats at inversion breakpoints should be considered on the lower end of the spectrum.

Characterization of gene content within inversions and their breakpoints

Based on the pan-genome inversion index, we identified a total of ~971 annotated genes per genome within inversions or at their breakpoints (Supplementary Data 12). To investigate the possible effects of inversions on gene expression, we interrogated a set of transcriptome datasets derived from three O. sativa subpopulations (i.e., XI-adm: MH63, XI−1A: ZS97 and GJ-temp: Nipponbare (i.e., dataset#2 – see Methods). Based on a comparison of transcript abundance levels between Nipponbare and MH63 and ZS97 (across four tissue types – root, panicle, young leaf, and mature leaf) we detected that 5–12 genes from MH63 and 4–11 genes from ZS97 within inversions, 2–7 genes from MH63 and 2–4 genes from ZS97 in inverted flanking regions, and 19–42 genes from MH63 and 9–30 genes from ZS97 that were located in non-inverted randomly resampled 20 Kb regions, were differentially expressed (DEG, fold change >2, p value <0.01) (Supplementary Fig. 7).

To investigate the effect of inversions on the transcription of genes located at inversion breakpoints – i.e., about 55 genes per genome (Supplementary Data 12), we interrogated both our baseline RNA-Seq datasets (dataset#1- see Methods) and dataset#2 for changes in transcript abundance. On average, transcript evidence for 28 of the 55 genes per genome could be detected in the tissues tested (Supplementary Data 12). Of these, transcript abundance of an average of 20 genes per genome did not change due to the presence of duplicated genes at both ends of their inversion breakpoints (Supplementary Data 12). An example of this observation is represented by two OsNAS genes (NAS1 and NAS2) located at the breakpoints of INV0300350 (~4.3 Kb) (Fig. 6a, b). The remaining ~8 genes per genome were single copy and were disrupted the inversion events, leading to the absence of transcript evidence (Supplementary Data 12). For example, transcripts of the Nipponbare Fbox gene (Os11g0532600) could be detected in the four tissues tested. However, the first exon of this gene was disrupted in MH63 by INV1101460, resulting in transcript ablation (Fig. 6c, d).

Fig. 6: Transcript abundance of genes located at inversion breakpoints.
figure 6

a Two copies of the OsNAS gene lie at the ends of an inversion in the MH63 (XI-adm) genome. This inversion disrupted the 5’ UTR regions of two OsNAS genes (NAS1 and NAS2); b OsNAS transcript abundance in root tissue; c The coding sequence (CDS) of a Fbox gene was disrupted by an inversion in the MH63 (XI-adm) genome; d Fbox transcript abundance was suppressed in all tissues tested in the MH63 (XI-adm) genome.

Recombination rate and genomic inversions

To evaluate the effect of inversions on recombination frequency, a previously published recombinant inbred line (RIL-10) population of 210 inbred lines43 derived from a cross between O. sativa cv. XI-adm: MH63 and XI−1A: ZS97 was investigated. We detected 78 inversions between MH63 and ZS97, totaling 3.58 Mb and 3.51 Mb in size, based on the MH63 and ZS97 genome assemblies, respectively (Supplementary Data 13). The recombination rate along each chromosome was assessed by comparing genetic and physical distances between neighboring bins. The average recombination rate for each chromosome ranged from 5.95 (chromosome 6) to 9.92 (chromosome 12) cM/Mb, and varied from 0 to 153.93 cM/Mb across the genome with an average of 6.98 cM/Mb (Supplementary Fig. 8a). The average recombination rate over the 78 inverted regions was 4.00 cM/Mb (0–23.26 cM/Mb), which is significantly lower (Student’s t test, p = 0.0002) than that observed genome-wide (Supplementary Fig. 8b). These results indicate that a marked suppression of genetic recombination is associated with inversions.

Effect of large inversions on population SNP variation

The occurrence of inversions can affect DNA polymorphism at the population level in several ways, including increased divergence in the inverted region and changes in linkage disequilibrium (LD) patterns40. The latter is particularly interesting as it can affect SNPs that are mapped to positions megabases apart, and can be a confounding factor in LD-based analyses. To determine whether large O. sativa inversions (>100 Kb) left a trace in patterns of LD along the IRGSP RefSeq, we used the 3K-RGP dataset to examine LD blocks near inverted regions (n = 88) (Supplementary Data 14). An inversion fixed in a population may lead to the disruption of LD blocks, in which some SNPs flanking the inversion on one side are in LD with SNPs on the distal part of the inversion, but not on the adjacent part (Fig. 7), due to the reversed order of SNPs inside the inverted region in samples that carry the inversion allele. By an LD block we mean only a set of SNPs in high LD (r2 > 0.8 in this analysis).

Fig. 7: Population linkage disequilibrium analysis of large inversions.
figure 7

A schematic diagram of linkage disequilibrium (LD) block disruption arising from the presence of an inversion, as shown in A and B. a Cartoon view of an inversion with breakpoints disrupting two LD blocks; b Expected features of the corresponding LD heat map; c Example of SNP blocks in high LD that are disrupted by an inversion; d The panel shows alignments, with the inversion marked by dotted lines. Small vertical lines above the horizontal axis mark the location of SNPs constituting a disrupted LD block. Orange and blue colors delineate two LD blocks that are contiguous in the of GJ-trop1 population, but appear as split when aligned to the IRGSP RefSeq (GJ-temp). Disruption of Azucena (GJ-trop1) haplotype blocks along the IRGSP RefSeq in the region of INV030410, as shown in e and f; e Genotype heat map of the GJ-trop1 subpopulation (samples in rows, SNPs in columns; light yellow: reference call, orange: heterozygous, brown: homozygous variant); f LD heat map of the same subpopulation. Dotted lines show the inversion region. Darker colors show larger r2. Note that the scaling of X-axis in the genotype heat map is not uniform, allotting half of X-axis space to the inverted region. Source data are provided as a Source Data file.

Next, we examined the entire 3K-RGP variation data set and searched for LD blocks that connect the flanking regions of inversions, having no SNPs in the proximal parts of each inversion. Such blocks (Fig. 7) were found in nearly all large inversions (126 out of 147 [85.7%], or 74 out of 88 [84%] inversion clusters) (Supplementary Data 14) with exceptions falling into two categories: i.e., inversions in regions of complex chromosomal rearrangements (Chr04:14.1–15 Mb, Chr11:6.6–6.85 Mb, Chr11:9.4–9.7 Mb), and three putative recent inversions (INV0401400, INV0500310, INV1000040), each of which were found in single genomes and may lack sufficient frequencies in a population to contain traces of recombination. Some of the disrupted LD blocks contained a particularly large number of SNPs and were seen as a distinctive checkered pattern on LD heatmaps (Fig. 7). This comparatively large number of SNPs along with low haplotype diversity, despite the presence of recombination, could be a consequence of selective pressure.

Read more here: Source link