Introduction

Transposable elements (TEs) are a major component of eukaryotic genomes, particularly in plants ¹. Their ability to move, creating mutations by insertion/excision, and amplify in genomes - thus generating new copies that provide opportunities for recombination - make them a major source of genetic variability ². For this reason TEs are considered a major driver of plant genome evolution, both in the wild or under human selection ³. TE insertions and other structural variants, however, have been frequently overlooked when using genome wide association studies (GWAS) to look for genetic variants linked to interesting phenotypes ⁴. The recent development of efficient tools to identify transposon insertion polymorphisms (TIPs) ⁵ has allowed the incorporation of TIPs in these analyses, and several GWAS reports have shown that TIPs uncover additional variability linked to phenotypic traits ^6–8. Interestingly, TIPs often explain more phenotypic variance compared to SNPs, and can be used for genomic prediction ⁹. The reason for this could be that, as compared with SNPs, TE insertions are more frequently the causal mutation.

A major fraction of the mutations linked to crop domestication and breeding are associated with changes of the expression of genes involved in signal transduction ^10,11. TEs can alter gene transcription, activating or repressing them, by different means. On the one hand, insertion of TEs within a promoter can interfere with transcription, and the silencing of TEs inserted close to genes can result in repression of gene expression ^12–14. Alternatively, TE insertions can also result in gene overexpression ^15,16, as TEs bring their own transcriptional promoters that can induce expression of neighboring genes. This is particularly clear for LTR-retrotransposons which contain promoters in their 5′ and 3′ LTRs, and may provide nearby genes with alternative promoters. In addition, TEs can also contain transcription factor binding sites (TFBS) that can alter the expression of host genes ¹⁷; it has been shown, for example, that MITEs have frequently amplified and redistributed TFBS in plant genomes ¹⁸.

Here we explore the impact of the movement of TEs during the recent evolution of rice on the variability of gene expression in this species. We took advantage of a recently published transcriptional analysis of 208 rice varieties grown under wet paddy and intermittent drought conditions ¹⁹ to perform an expression quantitative trait locus (eQTL) GWAS using TIPs as genetic markers. We show that TIPs are frequently associated with changes of expression of rice genes and that many TIPs altering the expression of regulatory genes have been positively selected in indica and japonica rice subspecies.

Figure 1

Results

TIPs are associated with gene expression variation in rice

We used a recently published dataset containing genome resequencing and transcription data for rice (Oryza sativa) varieties to perform TIP-eQTL mapping analyses. Rice is one of the most important food crops in the world, and varieties are generally found in two subspecies – O. sativa ssp. japonica and ssp indica. Other minor variety groups include aus varieties, which appear closely related to indica, and basmati rice varieties that appear to be a hybrid between ssp. japonica and aus. The data we analyzed consists of an ‘indica’ dataset for 126 O. sativa ssp. indica and some circum-aus accessions, and 84 ‘japonica’ accessions comprising of O. sativa ssp. japonica and some circum-basmati varieties, described in Groen et al. (2020) ¹⁹. We predicted TIPs using Popoolation TE2 ²⁰, and identified 45,050 TIPs. Using the Minghui MH63 short read data and assembled draft reference genome ²¹ we estimated the performance of the TIP genotyping on this dataset to have 97.1% sensitivity and 92.2% precision (Supplementary text 1). Most of the TIPs detected are related to DNA transposons (84%), especially MITEs belonging to the Stowaway (24%) and Tourist (20%) superfamilies; LTR retrotransposons make up a minor fraction (9 %).

The 45,050 TIPs are distributed along the 12 chromosomes and 81.3 % of them are less than 5 kb away from an annotated gene (IRGSP RAP-DB gene models). This results in up to 81.4% of rice genes having a TIP less than 5 Kb distant. In order to analyze the potential of TIPs to generate gene expression variability in rice, we used the 3’ mRNA-seq expression data from leaves of adult plants grown in normally-watered soil (“wet” condition) and under intermittent drought-stress (“drought stress” condition) ¹⁹. We selected data from 15,549 genes showing expression in more than 10% of the samples. Due to the strong population structure present in rice we separately analyzed data in the indica (126 varieties) and japonica (82 varieties) datasets. We found a total of 563 significant associations between TIPs and gene expression levels in the indica population and 356 in the japonica population (TIP-eQTLs with MAF > 3%, FDR adjusted p value < 0.05 in at least one replicate) (Figure 1A and 1B, Supplementary Table 2). These involve 477 and 317 genes, in indica and japonica respectively, representing 3.1% and 2% of the expressed genes. The TIP-eQTLs were related to the different TE types (orders ²²), and the proportion of TE orders found in the significant TIPs did not differ significantly from their representation across the genome (Figure 1C).

In order to estimate the relative contribution of TIPs and SNPs to gene expression changes, we also performed eQTL mapping analyses using 1 million of randomly chosen SNPs as markers ¹⁹. We obtained 4,913 SNP-eQTLs in indica and 3,308 SNP-eQTLs in japonica associated with SNPs located less than 5 kb from the gene; each SNP-eQTL represents the leading SNP-gene association. This corresponds to ~31.6% of the genes in indica and ~21.3 % in japonica. For each gene with a significant association, we compared the proportion of variance (R²) explained by the most significantly associated TIPs and/or SNPs. For most genes that have an associated TIP, either there is no associated SNP or the TIP explains more expression variance than the associated SNP (~62% of the indica TIP-eQTLs and ~73% of the japonica TIP-eQTLs) (Figure 1D). This demonstrates that TIPs uncover genetic associations with changes of expression that are not observed when using SNPs as molecular markers.

The SNP-eQTLs found in our study are associated to both positive (47%) and negative (53%) effects on gene expression levels, very close to an even distribution. In contrast, TIPs more frequently (59%) have a negative effect. The proportion varies among TE superfamilies, with TIR/MULEs being close to 50%, and LTR/Gypsy having the most frequently negative impact on gene expression (67%, Supplementary Table 1). This suggests that, as compared with SNP-eQTLs, TIP-eQTLs are more frequently the likely causal mutation responsible for the change in expression, and that the effect of the insertion is more generally negative, in line with recent analyses done in Capsella ²³, especially for long elements such as LTR-retrotransposons.

An analysis of the TIP and SNP-eQTL location with respect to the associated genes shows that while SNP-eQTLs are more evenly distributed across the genic region (Figure 1F), TIP-eQTLs are over-represented in the 1 kb region upstream of the gene (Figure 1E), which is the region that most frequently contains promoter elements regulating transcription. This suggests that although some TIP-eQTLs could be in linkage disequilibrium with the causal mutation for the change of expression, an important fraction of TIP-eQTLs could be the actual causal mutation inducing expression variation.

It has been proposed that an important fraction of gene expression variation is deleterious, and therefore, alleles associated with major expression effects should be maintained at lower frequencies ^24,25. An analysis of the effect size of the different TIP-eQTLs taking into account their frequency in the population shows that, indeed, low and high-frequency TIPs (corresponding to rare alleles) show the highest positive or negative effect on the expression of the associated gene (Figure 2 A-D, Supplementary Figure 1). This is consistent with what is described as the rare variant effect ²⁵ where rare (and likely deleterious) mutations in a population are associated with greater effects on gene expression in a population.

Figure 2

A comparison of TIP-eQTLs associated with changes of expression in the two different growth conditions (wet and drought stress) shows that there is an important overlap of eQTLs between the two environments (38.5 % for indica TIP-eQTLs and 33.2% for japonica TIP-eQTLs) (Figure 2E). Such overlap increases to 90% (average of indica and japonica, wet) and 71 % (average of indica and japonica, drought stress) when we compare TIP-eQTLs of one condition with a non-FDR-corrected TIP-eQTLs in the other environment (Figure 2G). This result suggests that the majority of TIP-eQTLs are associated with changes of expression in both growth conditions, and that only a small number of TIP-eQTLs (81 in indica and 78 in japonica) are associated with stress-specific changes of expression. Nevertheless, among the TIP-eQTLs associated with changes of gene expression in both wet and drought stress, there are TIP associated with genes known to be involved in drought tolerance (Supplementary Table 2). As an example, TIP_49046 is associated with a reduction of the expression of the gene synaptotagmin-5 (OsSYT-5). OsSYT-5 encodes a Ca²⁺ sensing protein with a C2 domain, is expressed in both stressed and non-stressed plants, and it has been recently shown that its silencing enhances drought tolerance in rice ²⁶. Other examples are TIP_23764 and TIP_52367, associated with an increased expression of the gene encoding the OsSAPK10 ABA-activated protein kinase and the S-type euonymus-related lectin gene OsEULS2 in respectively, two genes known to mediate drought stress in rice ^27,28

TIP-eQTLs are present at different population frequencies in indica and japonica

Although the overlap of TIPs associated with changes in expression under wet and drought stress conditions is very high, the overlap between indica and japonica TIP-eQTLs is low, and only 90 TIP-eQTLs out of the 563 from indica and 356 from japonica are significantly associated with variation of gene expression in both subspecies (Figure 2E). Even when we make the less-stringent comparison of the TIP-eQTLs of one subspecies with the non-FDR-corrected TIP-eQTLs of the other, as much as ~69% of indica and ~60% of japonica TIP-eQTLs appear as subspecies-specific (Figure 2F). This result may suggest differences in the transcriptional networks in indica and japonica or, alternatively, that different alleles showing expression level differences may be differentially represented in the two subspecies. An analysis of the frequencies in indica and japonica of all the TIP-eQTLs described here shows that approximately one-third of the TIPs are only found in one of the two subspecies. This suggests that TEs have been actively inserting (or have been excised or eliminated from the population) very recently during rice evolution, likely post-domestication and during the crop diversification phase, as has already been proposed for some rice TE families ^29,30. Here we show that some of these recent TIPs are associated with changes in gene expression and may therefore have phenotypic consequences. We found 44 TIP-eQTLs corresponding to TIPs present only in japonica varieties and 64 TIP-eQTL corresponding to TIPs present only in indica varieties. Some of these insertions are found at high frequencies in its corresponding subspecies, which suggest that they may have been under positive selection (Supplementary Table 3).

However, two thirds (66%) of all TIP-eQTLs found in this study correspond to TIPs present in both subspecies, which suggests that these insertions are relatively old and were already present in the ancestor of indica and japonica rice. Remarkably, in most cases they are found at very different frequencies in the two subspecies (Figure 2I). This suggests that these insertions were already present in the ancestor of indica and japonica rice and have been differentially retained in the two subspecies. Interestingly, an analysis of all non-eQTL-TIPs shows that an important fraction is found at the same frequency in both genomes (Figure 2H), suggesting that the preferential retention of TIP-eQTLs in indica or japonica may be a specific phenomenon potentially linked to their impact on the expression of the associated genes.

We examined whether TIP-eQTLs present in both indica and japonica varieties are also found in either Oryza rufipogon and Oryza nivara, the wild rice relatives believed to be the ancestors of domesticated rice ³¹. We looked for the presence of the TIPs identified in rice in a set of 72 accessions of O. rufigogon and 10 accessions of O. nivara (collectively called rufipogon/nivara from now on) (Supplementary Table 4). Up to 552 of the 829 TIP-eQTLs present in indica and/or japonica (66%), can be found in rufipogon/nivara varieties, confirming that an important fraction of TIP-eQTLs found in indica and japonica may be relatively old insertions already present in the wild ancestors of domesticated rice. Not surprisingly, in most cases the TIP frequencies in rufipogon/nivara are different from the frequencies of these insertions in indica and/or japonica (Supplementary Table 5); these may have arisen either due to selection or genetic drift accompanying the bottleneck associated with crop evolution.

In order to identify possible selection of TIPs in indica and japonica populations, we used the Population Branch Statistic (PBS) method ³² which examines strong differentiation in frequencies between populations. This approach consists in comparing the three pairwise FST values between indica, japonica and rufipogon/nivara to estimate the frequency change in TIPs that occurred since the divergence of the two rice subspecies from its wild ancestor. We calculated PBS values for a total of 11,698 TIPs present in the japonica, indica and rufipogon/nivara populations, which included 354 TIPs that are eQTLs in in indica and/or japonica populations. TIP-eQTLs (Figure 3A) showed much higher PBS values in comparison to those not identified as eQTLs (Figure 3B) (mean absolute PBS of 0.06 vs 0.02, respectively, p < 0.01, Wilcoxon test). Furthermore, TIPs with the most extreme PBS values (above 95 percentile) are enriched for TIP-eQTLs (Fisher’s Exact test odds ratio = 4.87, p < 0.01), suggesting that a high fraction of TIP eQTLs have been positively selected in indica or japonica. Note that some of the TIPs not identified as eQTLs in this particular dataset may actually be linked to variation of gene expression in organs, developmental stages or environmental conditions different from the ones analyzed here. Interestingly, up to 35% of the TIP-eQTLs identified here, associated with 119 genes, fell within the top 10% of absolute PBS values for all TIPs (Fig 3A), suggesting that TIPs associated with gene expression variation are likely under differential selection during rice evolution.

Figure 3

We observe that TIPs can show evidence of selection in japonica, indica or both subspecies. As examples, Figure 3B and C shows the representation of the PBS analysis, and the frequency in the populations, of two TIP-eQTL whose frequency has greatly increased in indica (TIP_53500), or japonica (TIP_72732), compared with the mean PBS for all TIPs.

Gene variants selected during rice evolution: some examples

Genes linked with TIP-eQTLs showing extreme indica or japonica PBS metrics are good candidates for genes underlying adaptation between different rice subspecies. Among the genes with TIP-eQTLs that have extreme PBS values in indica or japonica we can find several examples that are known to regulate plant architecture, plant and grain development or abiotic stress responses.

Some of the most extreme PBS values are for TIPs associated with changes of expression of genes involved in the signal transduction of hormones, including brassinosteroids, ABA, ethylene, jasmonic acid (JA) and auxin (Supplementary Table 5). We found four different TIP-eQTLs with high PBS values associated with the EG2/OsJAZ1 gene (Os04g0653000), a locus encoding a JA signaling repressor ³³; the four insertions are related to MITEs of the Tourist transposon superfamily. Three of the insertions (TIP_32891, TIP_32892 and TIP_32894), located upstream (3.4 kb and 2.7 kb) and downstream (150 bp) of the gene, are associated with an increase of expression of EG2, and are physically linked to one another (mean r² value = 0.88). In contrast, the fourth insertion (TIP_32893), is not linked to them, is within the first intron of the gene and is associated with expression reduction.

The recent advances in characterizing the pangenome of rice and the super-pangenome that includes its wild ancestors ^34–36 has allowed us to characterize the locus in indica, japonica and rufipogon/nivara. This analysis showed that EG2 is primarily found in two different haplotypes (Figure 4A). One haplotype, Hap A, contains the three MITE insertions associated with higher expression of EG2 (Figure 4A, 4B, and Supplementary Figure 2A), while the second haplotype, Hap B, contains the fourth MITE insertion and is associated with reduced EG2 expression (Figure 4A, 4B and Supplementary Figure 2A). Hap A was identified to be at high frequency in rufipogon/nivara (48%), whereas Hap B is at lower frequency (16%). The frequency of both Hap A and Hap B is slightly increased in japonica (65% and 28%, respectively). In indica, however, there is a clear reduction in the frequency of Hap A and a concomitant strong frequency increase of Hap B, associated with reduced expression in indica (16% and 83%, respectively). Interestingly, EG2/OsJAZ1 is a repressor of spikelet development ³³ and the number of differentiated spikelets per panicle tends to be lower in japonica as compared with indica ³⁷.

Figure 4

We looked more closely for signs of selection at this specific locus, by examining levels of nucleotide diversity (π³⁸) and diversity in haplotypes homozygosity tracts (H12 ³⁹). Consistent with the expression (and possible spikelet phenotypic differences), we found that in japonica the haplotype associated with increased EG2 expression (Hap A) was under positive selection, whereas the haplotype associated with decreased expression of EG2 (Hap B) did not show any sign of selection; this indicates that positive selection acts on individuals carrying the haplotype for high expression of EG2 in japonica (Figure 4C-D). Specifically, we observe significantly lower π and higher H12 in Hap A, associated with increased EG2 expression, while this pattern was not significant for Hap B (supplementary Table 6). However, in indica the opposite pattern appears to prevail, wherein there is a stronger evidence of selection for Hap B, but less so for Hap A (Figure 4E-F). This suggests positive selection towards higher expression of EG2 in japonica, and lower expression of EG2 in indica, which could potentially explain the higher number of spikelets per panicle observed in the latter³⁷.

Another good candidate for genes underlying adaptation of rice is the OsGAP gene. There are two different TIP-eQTLs with high PBS values associated with changes of expression of the OsGAP gene (Os07g0500300) in japonica, whereas they do not correlate with changes in expression in indica (Supplementary Figure 2B). OsGAP encodes a putative GTPase activating protein, similar to CAR proteins (C2-domain abscisic acid-related proteins), that play an important role in ABA signal transduction in Arabidopsis ⁴⁰. The first TIP (TIP_50057) is located ~4 kb upstream of the OsGAP gene and is associated with an increase of expression of the gene in japonica whereas the second insertion (TIP_50059) is located within the first intron and is associated with a decrease in expression in japonica (Figure 5A-B). The analysis of the locus using the super-pangenome data shows that OsGAP is found in two different haplotypes, one containing TIP_50057 (Hap A) and the other containing TIP_50059, together with some additional structural differences (Hap B; Figure 5A). The two TIPs defining the two haplotypes are present in the rufipogon/nivara population at complementary frequencies (73% and 19%, respectively). Interestingly, the proportion of the two haplotypes in cultivated rice seem to be very different, with an increased frequency of the haplotype associated with a reduced expression of OsGAP (Hap B) in both japonica (~50%) and indica (~95%) were it reaches near fixation.

Figure 5

It has been proposed that OsGAP is a negative regulator of ABA signaling in seed germination and dormancy, and that reduced expression of OsGAP may prevent rice pre-harvest sprouting (PHS) ⁴⁰. Reduced seed dormancy is a common target for selection of cultivated varieties, but this trait may come at a cost of a higher PHS risk, which may be a problem, especially in regions where heavy rain is common during the harvest season. Therefore, the appropriate degree of dormancy may depend on the agroecological conditions in which a particular variety is grown. We see here that two independent TIPs are associated with changes of OsGAP expression in japonica but not in indica (Figure 5A, Supplementary Figure 2B). Interestingly, there are clear and strong signs of selection for increased expression of OsGAP in japonica, wherein Hap A was under positive selection, and Hap B was selectively neutral, indicating positive selection acts on only those individuals carrying haplotype for high expression of OsGAP in japonica (Figure 5C,D; Supplementary Table 6). In contrast, there is no evidence for selection in indica varieties despite the near fixation of one haplotype. This suggests that an increased level of OsGAP may be relevant for the control of dormancy in japonica but not in indica (despite its near-fixation in the latter), which would be in line with recent data suggesting that seed dormancy is regulated by different genes/alleles in indica and japonica ⁴¹.

Finally, we find a TIP-eQTL whose frequency has greatly increased in cultivated rice with respect to rufipogon/nivara. TIP_45706, corresponding to an insertion of a gypsy-like LTR-RT at ~1 kb upstream of the OsMPH1 gene (Figure 6A). This insertion is associated with a decrease in the expression of OsMPH1 in indica and possibly japonica varieties (Figure 6B, Supplementary Figure 2C), and is found at low frequency in rufipogon/nivara (~11%) but at high frequency in both indica (70%) and japonica (95%). OsMPH1 encodes a MYB-like transcription factor that has been shown to regulate plant height, its reduced expression resulting in shorter plants ⁴². A reduction in plant height due to a mutation in the SD1 gene, which encodes the gibberellin biosynthesis gene GA-20ox, was at the origin of the Green Revolution in the 1960s, but it has been shown that alleles of SD1 resulting in shorter culm length were also selected during the domestication of japonica rice ⁴³. Our results suggest possible parallel selection of alleles in other genes that may lead to similar phenotypes. This is further reinforced by strong selection acting on the haplotype containing the TIP_45706 insertion in both indica and japonica (Figure 6C and 6D Supplementary table 6), leading to reduced expression of OsMPH1 and thus shorter plants.

Figure 6

Discussion

Transposons are a major source of genome variability, and are known to affect gene expression in numerous ways. The importance of rare variants and unfixed TE insertions for gene expression variation in humans ⁴⁴ and plants ²³ has already been demonstrated. A recent pangenome analysis in tomato has shown that the phenotypic variation of important crop traits is linked to structural variants present in the species, most of which are actually TIPs, which are associated to subtle changes of transcription of genes involved in signal transduction ⁴⁵. In the present study we evaluated the importance of TIP-related expression variability in the recent evolution of rice. To this end we performed a TIP-eQTL mapping using expression data from rice varieties from the O. sativa ssp. indica and O. sativa ssp. japonica subspecies. We show here that using TIPs in addition to SNPs as genetic information allows the uncovering of additional genomic associations to changes in gene expression, and that when both TIPs and SNPs are both associated, TIPs often explain more of the variance in expression. This is in line with recent reports that have also used TIPs for GWAS with different phenotypes in tomato and rice ^7,8, and suggests that TIPs may be more frequently the causal mutation of a particular phenotype than SNPs. The concentration of TIP-eQTLs in upstream regions of genes, where most transcriptional regulatory elements usually reside, also suggest that a high fraction of the TIP-eQTLs here described are the actual mutation underlying gene expression variation.

The close association of some TEs with genes could also partially explain the frequent association of more than one TIP with changes of expression of particular genes. Indeed, among the 718 genes with TIP-eQTLs described here, 18% have two of more TIPs associated with changes of expression. However, not all genes have the same number of TIPs. Among the 30 genes with TIPs in indica and japonica that have the highest population branch statistics (PBS) – an indicator of strong differential selection - 13 have more than one TIP-eQTL. This suggests that TIPs are a major source of gene expression variability in rice, in particular in genes that may have been strongly selected during the recent evolution of rice. In some cases, the different TIPs linked to a gene are associated with opposite effects on its expression and are present in different haplotypes. In these instances, as shown for the OsJAZ1 and OsGAP genes, the two haplotypes have strong signs of selection, although one positively and the other negatively, pointing to role these TE insertions may play in adaptive evolution of these crop subspecies.

The results presented here show that most TIPs associated with variation of gene expression in indica and/or japonica are relatively old insertions that are also present in rufipogon/nivara, the wild ancestors of rice. This is in line with recent data showing that a high number of the structural variants associated with changes of expression in tomato were already present in its wild ancestor ⁴⁵ and highlights the importance of the standing variation, already present in the wild ancestors, for crop adaptation. It is assumed that most TE insertions are selectively neutral or slightly deleterious ⁴⁶, and the presence of TIPs associated with expression variation of genes in the wild ancestors of rice, suggests that this variation can be well tolerated in wild rice. Interestingly, we show here that many of these TIPs have been positive or negative selected in rice populations, which shows that they translate into selectable phenotypic differences in the agroecological conditions of cultivated rice. Indeed, we show examples of selected variants with modified expression of genes known to be linked to important traits that were targets of selection. It has been already proposed that a major fraction of the mutations linked to crop domestication and breeding are associated with changes of gene expression involved in signal transduction ^10,11. Here we show that specific expression variants of genes involved in signal transduction have been differentially selected in indica and japonica rice populations. In addition, our results also point to TEs as a major driver of gene expression variation selected during crop adaptation and breeding.

Methods

TIP detection

Re-sequencing data for 126 indica and 82 japonica rice accessions was obtained from Groen et al ¹⁹ (Bioproject accessions PRJNA557122, PRJNA422249 and PRJEB6180). BBDuK (https://sourceforge.net/projects/bbmap/) was used for adaptor and quality trimming. Clean reads were aligned to the Nipponbare reference genome ⁴⁷ using BWA ⁴⁸. PoPoolationTE2 ²⁰ was used to detect TIPs in the mode “joint” using Nipponbare TE annotation described by Ou et al ⁴⁷. TIPs with a zygosity lower than 0.25 in all samples were excluded to avoid false positives. TIPs were further filtered using the parameters –-min-count 5, --max-otherte-count 2, --max-structvar-count 2, and only those having MAF higher than 3% and no missing data were kept. Finally, the TIP matrix was transformed to binary form using zygosity cutoff of 0.05 to define an insertion as present.

TIP and SNP-eQTL mapping

Transcriptome data for the 208 rice accessions (wet and drought stress conditions; tree independent replicates per condition) was obtained from Groen el al ¹⁹ and separated into individual replicate matrices. Transcripts expressed in more than 10% of the samples on each replicate were extracted and counts were normalized using the yst function of DESeq2 ⁴⁹. TIP and SNP-eQTL mapping were performed using Matrix eQTL software ⁵⁰, applying the simple linear regression model and including subpopulation groups as covariates. We used a cut-off distance of 5,000bp to identify cis-eQTLs and a 5% False Discovery Rate threshold for multiple testing corrections.

Selection analyses

For the Population Branch Statistics (PBS) analysis we looked for the presence of rice TIPs in the 82 accessions belonging to the rufipogon/nivara population (Supplementary Table 4) and retained only those present in the three populations (indica, japonica, rufipogon/nivara), resulting in a matrix of 13,622 TIPs and 288 accessions. The TIP matrix was filtered to remove TIPs with more than 5% missing data, or with MAF < 5%. The remaining missing data (1.7%) was imputed using the “wright” algorithm implemented in SNPready. The clean matrix (11,698) was used to calculate the individual PBS values per TIP, following the formulas described in Yi et al., (2010)³².

Along with PBS, to evaluate whether TIP-eQTLs show signs of selection, we estimated nucleotide diversity (pi ³⁸), which is a SFS (site-frequency spectrum) measure to test for the presence of selection. This was done using a custom script in a window of 25SNPs with a sliding window of 5. Further, we also estimated the H12 homozygosity statistic ³⁹ that can detect the presence of both hard and soft sweeps associated with selection. Since we expect long runs of homozygosity in the regions around the selected loci, H12 was calcuted in windows of 50 SNPs, with a sliding window of 25. Both these statistics were estimated separately for indica and japonica, and for the samples with TIP insertions present and absent. This was for a range of 100kb (50kb up-and downstream the TIP insertion), using SNPs identified from Groen et al ¹⁹ but without the 1000bp linkage thinning. To test for the significance of the statistics, we performed a permutation test (N=1000) with random sampling with the same number of individuals that were used to estimate the selection statistics. The statistics were estimated for these 1000 random pulls in the same regions, and using the same window sizes as the original statistics.

Acknowledgements

We thank Jae Yong Choi, Simon C. Groen and Sebastián Ramos-Onsins for helpful discussions. Noemia Morales-Diaz holds Predoctoral fellowship PRE2020-095111-funded by MCIN/AEI/ 10.13039/501100011033 and by “ESF Investing in your future”. Raúl Castanera holds a Juan de la Cierva Incorporación Postdoctoral fellowship IJC2020– 045949-I funded by MCIN/AEI/10.13039/501100011033 and by the “European Union NextGenerationEU/PRTR”. This work was funded by the Spanish Ministerio de Ciencia e Innovación (PID2019-106374RB-I00/AEI/10.13039/501100011033). We also acknowledge financial support from the Spanish Ministerio de Economía y Competitividad through the “Severo Ochoa Program for Centers of Excellence in R&D” 2016–2019 (SEV-2015-0533) and CEX2019–000902-S funded by MCIN/AEI/10.13039/501100011033.

This work was also supported by grants SEV - 2015 - 0533 and CEX2019-000902-S funded by MCIN/AEI/10.13039/501100011033, and by the CERCA Program / Generalitat de Catalunya.

Author contributions

RC and JMC conceived the project. RC and JMC planned the experiments, together with SG and MDP for the selection analyses. RC performed most experiments with the help of NMD, who also solved the structure of the different haplotypes. SG performed the analyses of selection. All authors discussed the results JMC and RC drafted the article. All authors contributed to the paper.