Leveraging ancestral recombination graphs for quantitative genetic analysis of rice yield in indica and japonica subspecies

Ancestral Recombination Graph; Genomic prediction; Rice16 Introduction1 R ice (Oryza sativa L.) is a major food staple, with two main2 subspecies, indica and japonica. Despite their co-existence3 in many regions, indica and japonica are generally treated as4 separate populations in breeding programs. Combining these5 populations could improve the efficacy of current genomic pre-6 diction and genome-wide association studies; however, this7 requires appropriate modelling of the populations’ genetic his-8 tory . Ancestral recombination graph (ARG) is the ultimate way9 to represent the complete genetic history of individuals and10 their genomes. ARGs have been shown to be effective in many11 population genetic studies, and hold great potential to improve12 quantitative genetic modelling. In this study , we evaluate the use13 of an ARG, encoded as tree sequences, to improve quantitative14 genetic modelling in indica and japonica rice subspecies.15 Genetic analysis of cultivated rice varieties revealed two main16 subspecies, indica and japonica (Glaszmann 1987; Garris et al.17 2005), which share a common ancestor with the wild O. rufi-18 pogon species (Huang et al. 2012). They were estimated to have19 diverged ∼0.4 million years ago (Zhu and Ge 2005), by a single20 domestication event with subsequent divergence and migra-21 tion between subspecies, or by multiple domestication events22 (Caicedo et al. 2007; Liu et al. 2024). Due to this deep divergence, 23 the cultivation, breeding, and quantitative genetic analyses are 24 generally conducted separately for indica and japonica. However, 25 in some countries in Asia and elsewhere, such as Uruguay in 26 South America, indica and japonica coexist (Martínez et al. 2014). 27 In these places, small- and medium-scale rice breeding programs 28 of indica and japonica occur, and individual populations of each 29 subspecies have limited size. This raises the possibility of consid- 30 ering joint quantitative genetic analyses and breeding operations 31 for both subspecies, but the genetic history of rice subspecies 32 requires appropriate modelling (Caicedo et al. 2007; Liu et al. 33 2024). 34 Rice breeders generally select and improve their populations 35 based on estimated genetic values, which are obtained from a 36 linear mixed model with a relationship matrix (RM) describing 37 the genetic similarity between individuals (Henderson 1984). 38 The RM can be built from different sources, such as pedigree 39 or genomic data. A pedigree-based relationship matrix (PRM) 40 expresses the expected identity by descent between pairs of in- 41 dividuals relative to the base population of a known pedigree 42 (Henderson 1984; Mrode 2014). In contrast, the genomic RMs 43 proposed by VanRaden (2008) express the realized identity by 44 state across single nucleotide polymorphism (SNP) sites. This 45 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 20, 2025. ; https://doi.org/10.1101/2025.01.14.633033doi: bioRxiv preprint 2 Genomic Prediction with Ancestral Recombination Graphs in rice notion captures the relationships between pairs of individuals1 due to shared mutations that segregate between and within2 families. We refer to this type of RM as a site-based relation-3 ship matrix (SRM), because it is based on current variation at4 polymorphic sites typed by SNP markers. The SRM enhances5 quantitative genetic modelling in terms of genomic prediction6 relative to the PRM and enables genome-wide association stud-7 ies (de los Campos et al. 2013; Yu et al. 2006). However, the8 commonly used SRM does not fully reflect the genetic history of9 sampled individuals because it does not capture all SNP sites,10 and origin of mutations (Young 2022; Fan et al. 2022). This can11 reduce the efficacy of these applications for multiple (divergent)12 populations, like indica and japonica rice, so there remains the13 need to consider appropriate RMs for joint quantitative genetic14 analyses of these subspecies.15 Joint quantitative genetic analysis of multiple populations16 is an active area of research (Martin et al. 2017; Ding et al. 2023;17 Mester et al. 2023; Warburton et al. 2023; Zhang et al. 2023b; Sun18 et al. 2024). One of the key considerations for joint analyses of19 small populations is that the construction of a sufficiently large20 training dataset is likely to include diverged populations with21 specific recombination and mutation events, which can coun-22 teract the advantages of constructing a larger training dataset23 (Lorenz and Smith 2015). Therefore, the joint quantitative genetic24 analysis of indica and japonica could be beneficial, but requires a25 better representation and modelling of the quantitative genetic26 diversity between and within the subspecies.27 An ARG describes the genetic history of sampled genomes28 by encoding their shared ancestry , mutation, and recombina-29 tion events (Griffiths and Marjoram 1996; Brandt et al. 2024;30 Lewanski et al. 2024; Nielsen et al. 2024; Wong et al. 2024). ARGs31 have been used to recover relationships between individuals32 and populations (e.g., Kelleher et al. 2019; Wohns et al. 2022; Fan33 et al. 2022), estimate coalescence times (e.g., Brandt et al. 2022),34 map quantitative trait loci (QTL) (e.g., Link et al. 2023), and esti-35 mate heritability and genome-wide association (e.g., Zhang et al.36 2023a). Within an ARG, shared ancestry at a locus is encoded37 by a local tree, where the nodes represent the sampled haploid38 genomes and their ancestors (Brandt et al. 2024; Lewanski et al.39 2024; Nielsen et al. 2024; Wong et al. 2024). The nodes are con-40 nected with edges, represented as branches in a local tree, which41 indicate the lineages of descent between haploid genomes. Mu-42 tations on the branches generate variation between genomes at43 a given locus. Recombinations between loci change the lines of44 descent between haploid genomes and hence change the topol-45 ogy of local trees along the genome (Brandt et al. 2024; Lewanski46 et al. 2024; Nielsen et al. 2024; Wong et al. 2024). Recent advances47 in methodology and software implementations now enable the48 inference of ARGs for many genomes and populations (Speidel49 et al. 2019; Kelleher et al. 2019; Zhang et al. 2023a; Harris 2023;50 Gunnarsson et al. 2024; Nielsen et al. 2024)51 In this study , we use the tree sequence representation of an52 ARG to improve quantitative genetic modelling of indica and53 japonica rice subspecies. A tree sequence succinctly encodes a54 sequence of local trees along the genome (Kelleher et al. 2016,55 2018). Tree sequences enable efficient genome storage in simula-56 tions (Kelleher et al. 2016, 2018; Haller et al. 2019) and population57 genetic studies (Kelleher et al. 2019), as well as fast calculations58 of population genetic statistics (Ralph et al. 2020). There is now59 a collection of software to work with tree sequences, such as60 inferring tree sequence from phased genomes (Kelleher et al.61 2019), dating nodes of a tree sequence (Wohns et al. 2022), and62 calculating population genetic statistics from a tree sequence 63 (Ralph et al. 2020). Importantly , the work of Ralph (2019) and 64 Ralph et al. (2020) provide a framework for calculating popu- 65 lation genetic statistics based on site and branch information, 66 including a site-based and branch-based RM from an ARG (B. 67 Lehmann et al., personal communication, November 27, 2020). 68 We used tskit implementation of branch-based RM (Kelleher et al. 69 2021a), which we call BRM. The BRM captures the rich genetic 70 history of sampled individuals, including untyped sites, link- 71 age between sites, and enables estimation of haplotype effects. 72 Recently ,Fan et al. (2022), Zhang et al. (2023a), Link et al. (2023), 73 Gunnarsson et al. (2024), and Zhu et al. (2024) showed that the 74 BRM captures the genetic diversity better than the SRM, opening 75 opportunities to improve quantitative genetic modelling across 76 multiple populations, such as rice subspecies. This work aimed 77 to evaluate the use of tree sequences as a tool to build a RM for 78 a joint quantitative genetic analysis of indica and japonica rice 79 genotypes from a breeding program in Uruguay . To this end, we 80 inferred and analyzed tree sequences for a sample of individu- 81 als from the Uruguayan National Rice Breeding Program, built 82 the BRM, and used it for genomic prediction and genome-wide 83 association study of grain yield. The results of this study can 84 guide strategies for joint selective breeding and management 85 of genetic diversity in indica and japonica rice and other species 86 with a similar population structure. 87

88 We used genomic, pedigree, and phenotypic data from the 89 Uruguayan National Rice Breeding Program of indica and japon- 90 ica rice between 1997 and 2020 (Martínez et al. 2014; Rebollo 91 et al. 2023a). We first inferred the rice ancestral alleles and the 92 ARG for all individuals using a tree sequence with dated nodes. 93 We then analyzed the properties of inferred and dated tree se- 94 quences, including genealogical nearest neighbor analysis, and 95 inspected the local trees at two selected regions. Using differ- 96 ent cross-validation (CV) scenarios, we evaluated the ability of 97 linear mixed models using either the PRM, SRM, or BRM to 98 predict grain yield between and within the subspecies. Finally , 99 we inspected estimated SNP site effects to understand how the 100 estimates map onto branches of the ARG. 101 Data 102 Genomic, pedigree, and phenotypic data were retrieved from 103 the Uruguayan National Rice Breeding Program. A detailed 104 description of the data and their availability is given in Rebollo 105 et al. (2023a), and summarised in the following. 106 Genomic data were available for 965 late-stage inbred lines of 107 indica (395 individuals) and japonica (570 individuals), generated 108 with Genotyping-by-Sequencing (GBS) of ApeKI enzyme-cut 109 DNA fragments (Elshire et al. 2011; Rebollo et al. 2023a). We 110 aligned the sequences to the Nipponbare reference genome MSU 111 version 7.0 (Kawahara et al. 2013) using bwa (Li and Durbin 112 2009). Next, we used the TASSEL 5.0 pipeline (Bradbury et al. 113 2007) to call SNP sites and individual genotypes. We performed 114 quality control and retained SNP sites with minor allele fre- 115 quency ≥ 0.03, missing data ≤ 50%, and observed heterozygos- 116 ity ≤ 15%. The resulting genotypes were phased and imputed 117 with BEAGLE version 5.1 ( Browning et al. 2018). We performed 118 all steps above for both subspecies together (ALL), and sepa- 119 rately for indica (IND) and japonica (JAP). A total of 61,260 SNP 120 sites were retained for ALL, 50,854 for IND, and 23,614 for JAP . 121 Some SNP sites were shared between genomic datasets, while 122 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 20, 2025. ; https://doi.org/10.1101/2025.01.14.633033doi: bioRxiv preprint Rebollo et al. 2025 3 some were exclusive to each dataset (Supplementary Figure 1).1 In line with the aims of the study , the majority of the analyses2 and subsequent results focus on the ALL genomic dataset.3 Pedigree data were available for 915 of the 965 individuals4 with genomic data (95%), producing a total of 1,207 individuals5 in the pedigree including 292 ancestors. There were up to six6 generations of ancestors with 204 founder individuals. Individu-7 als without pedigree data were included as nominally unrelated8 to the remaining individuals, i.e., treated as founders. The num-9 ber of self-pollination generations were stored with the pedigree10 for later calculation of the PRM.11 Grain yield data were available for 936 of the 965 individuals12 with genomic data (97%), which are a subset of the 19,447 indi-13 viduals tested across 23 years of field evaluation (Rebollo et al.14 2023a). The selected grain yield data came from 846 field trials15 conducted between 1997 and 2020 in three locations: Paso de16 la Laguna (33.27 S, 54.17 W) with 743 trials; Paso Farías (30.5417 S, 57.26 W) with 92 trials; and Pueblo del Barro (31.93 S, 55.3818 W) with 11 trials. Each trial had a randomized complete block19 design with two to four blocks. All trials were conducted under20 current production management standards for irrigated condi-21 tions (Batello et al. 2013). Grain yield was recorded from 1.20 m22 x 2.0 m plots (2.4 m 2) and converted to kg/ha. The total number23 of recorded plots was 23,311. The number of individuals present24 each year ranged from 12 to 653 (mean of 173), while the number25 of individuals shared between pairs of years ranged from 3 to26 640. Field trials were grouped into 49 environments, defined27 by their year location combination. Each environment had 128 to 48 field trials (mean of 17), and 10 to 653 individuals with29 genomic data (mean of 95). Environments were excluded if they30 contained less than 20 individuals. The final genomic, pedigree,31 and phenotypic datasets included 22,741 records across 40 envi-32 ronments for 936 individuals, with 381 indica and 555 japonica33 individuals.34 Ancestral recombination graph35 We combined our genomic dataset and Ensembl’s alignment36 of Oryza sativa with other Oryza spp. (Yates et al. 2021) to infer37 the rice ancestral alleles using O. rufipogon, O. glaberrima, and38 O. meridionalis as external species. We converted Ensembl Mul-39 tiple Format (EMF) to Multiple Alignment Format (MAF) and40 accessed specific regions in the external species corresponding to41 SNP sites in our genomic dataset using WGAbed (Corcoran and42 Barton 2021). We inferred the rice ancestral alleles using est-sfs43 (Keightley and Jackson 2018) with Rate-6 model of nucleotide44 substitution and 10 random starting values. We declared the45 most probable allele at the common ancestor with the first out-46 group (O. rufipogon) as the ancestral allele. Where we could not47 infer the ancestral allele, we declared the major allele as ances-48 tral. We summarised this information by computing the total49 number of SNP sites inferred, and the number and percentage50 of times that the ancestral allele was the major, minor, or another51 allele. Finally , we computed and compared the alternative (i.e.,52 non-ancestral) and minor allele frequency spectrum of the ALL,53 IND, and JAP genomic datasets.54 We inferred an ARG encoded as a tree sequence for each chro-55 mosome with tsinfer 0.2.1 (Kelleher et al. 2019, 2021b). The56 succinct encoding of tree sequences is achieved by a collec-57 tion of core tables that describe: (i) haploid genomes (nodes)58 and their date of existence, (ii) shared DNA between nodes59 (edges/branches), (iii) SNP sites with their ancestral state, and60 (iv) mutation events on the edges giving rise to polymorphisms61 at the sites, alleles of the nodes, and genotypes of individuals. 62 For inferring tree sequences, we used a recombination rate of 63 4.53 × 10−8 per base pair per generation (Si et al. 2015), the de- 64 fault mismatch ratio parameter of 1, and physical SNP positions 65 from the Variant Call Format (VCF) file. Finally , we summarised 66 information on the number of nodes, edges, trees, SNP sites, and 67 mutations for the tree sequence of each chromosome. 68 We dated the ARG by inferring the age of the nodes with 69 tsdate 0.1.4 (Wohns et al. 2022). For this inference, we used a 70 mutation rate of 2.2 × 10−9 per base pair per generation (Yang 71 et al. 2015). This dating required knowledge of the effective 72 population size (Ne) and as there are contrasting estimates of Ne 73 depending on the type of data and time, we used three different 74 values (23, 1,500, and 150,000) to evaluate its impact on dating 75 and downstream analyses. We estimated the current Ne of 23 76 from the pedigree-based rate of coancestry (Pérez-Enciso 1995), 77 which is concordant with Rutkoski (2019). We estimated the 78 recent Ne of 1500 from the linkage disequilibrium of the ALL 79 genomic dataset at 50 generations ago with GONE 29.08.2021 80 (Santiago et al. 2020). We ran GONE with default parameters, 81 except for the recombination rate mentioned above and the max- 82 imum number of recombination bins set to 0.005. We used the 83 ancient Ne of 150,000 from the work of Caicedo et al. (2007). Re- 84 sults presented in the main text correspond to Ne = 1500 and we 85 show the sensitivity of the results to the other two Ne values in 86 Supplementary Material. From the dated ARG, we estimated 87 the age of SNP sites as the age of the first ancestor above which 88 we inferred a mutation. Finally , we summarized the information 89 on the age of ancestors, mutations, and SNP sites. 90 To extract biological signals from the ARG, we calculated the 91 genealogical nearest neighbors (GNN) statistic (Kelleher et al. 92 2019) for all individuals and inspected two local trees of interest. 93 The GNN describes the topology of an ARG with respect to 94 a reference set (subspecies), to describe the identity of nearest 95 neighbors within a local tree and summarized across the local 96 trees (Kelleher et al. 2019). We calculated a GNN matrix for 97 each chromosome between all pairs of individuals with tskit 98 0.3.7. We then combined these matrices into a single GNN ma- 99 trix by summing across chromosomes and weighting by their 100 relative proportions of the genome. We hierarchically clustered 101 the individuals with the average method using SciPy (Virtanen 102 et al. 2020). The first local tree of interest was between positions 103 32,645,584 and 32,648,666 bp of chromosome three, covering the 104 drought and salt tolerance (DST ) gene, which is associated with 105 panicle length in japonica only (Bai et al. 2016). The second local 106 tree of interest was between positions 31,583,546 and 31,596,967 107 bp of chromosome four, covering position 31,590,530 bp of chro- 108 mosome four which is associated with the number of panicle 109 secondary branches (NSB) and number of secondary spikelets 110 per secondary branch (NSSB) in both indica and japonica (Bai et al. 111 2016). 112 Relationship matrices 113 We calculated the PRM considering the number of generations 114 of self-pollination with preGSf90 (Aguilar et al. 2014; Rebollo 115 et al. 2020). The SRM was calculated with the rrBLUP R pack- 116 age (Endelman 2011) as MM⊤/c, where M is the SNP genotype 117 matrix centered by the mean allele dosage at each site and c 118 is the sum of the expected heterozygosities across SNP sites 119 (VanRaden 2008). A BRM was calculated for each chromosome 120 from the inferred tree sequence by computing the total area of 121 shared branches between and within each pair of homologous 122 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 20, 2025. ; https://doi.org/10.1101/2025.01.14.633033doi: bioRxiv preprint 4 Genomic Prediction with Ancestral Recombination Graphs in rice chromosomes of all individuals, relative to the remaining ho-1 mologous chromosome pairs using tskit 0.3.7 (Kelleher et al.2 2021a) (B. Lehmann et al., personal communication, November3 27, 2020). These matrices were then combined into a single BRM4 by summing across chromosomes and weighting by their rela-5 tive proportions of the genome. Branch-based statistics, such as6 the BRM, are calculated from the length and span of branches,7 giving the expected value of its dual site-based statistic under8 the infinite-sites mutation model (Ralph 2019; Ralph et al. 2020).9 For the purpose of subsequent quantitative genetic analyses, we10 write the BRM as MDM⊤, where M is the SNP genotype matrix11 and D is a relationship matrix between SNP site effects, induced12 through the tree sequence structure (see the Appendix).13 Quantitative genetic analysis14 We fitted various linear mixed models, with a particular focus15 on modelling different structures for genotype-by-environment16 (GxE) interaction and different genotype relationship matrices.17 The phenotypic dataset comprised v = 936 individuals evalu-18 ated in t = 828 trials across p = 40 environments (year locations)19 with n = 22, 741 records in total. Let the n-vector of phenotypic20 records be given by y = ( y⊤ 1, . . . , y⊤ p)⊤, where yj is the nj-vector21 for the jth environment. The linear mixed model for y is given22 by:23 y = Xb + Zuu + Zpp + e, (1) where b is a vector of fixed effects with design matrix X, u is the24 vp-vector of random genetic values, ordered as individuals in25 environments, with design matrix Zu, p is a vector of random26 non-genetic effects with design matrix Zp, and e is the n-vector27 of residuals. The fixed effects included the overall and environ-28 ment means, while the random non-genetic effects captured the29 blocking structure of trials within environments. The random30 effects and residuals were assumed to be mutually independent31 following a multivariate normal distribution.32 We fitted seven linear mixed models, as summarised in Sup-33 plementary Table 1. All models included diagonal covariance34 structures with separate block and residual variances for each35 environment. Note that the absence of field layouts precluded36 fitting residual spatial models. All models also included a37 separable covariance structure for the genetic values given by38 var(u) = Ge ⊗ G, where Ge is an unknown covariance structure39 between environments and G is a known RM between individ-40 uals. We considered seven forms of Ge to model different GxE41 interaction patterns and three forms of G to model the different42 sources of genetic data (i.e., the PRM, SRM, and BRM). Details43 on the different forms of Ge can be found in Tolhurst et al. (2022),44 but briefly:45 • Model 1 assumes independent genetic values between en-46 vironments, modelled by a diagonal covariance structure47 with a separate genetic variance for each environment.48 • Model 2 assumes correlated genetic values between envi-49 ronments, modelled by a uniform covariance structure with50 a single genetic variance and covariance.51 • Model 3 is an extension of model 2 that includes an addi-52 tional interaction variance across environments, producing53 a compound symmetric covariance structure.54 • Model 4 is an extension of model 3 that includes a separate55 interaction variance for each environment, rather than a56 single variance across environments.57 • Models 5, 6, and 7 include factor analytic covariance struc-58 tures of order 1, 2, and 3, respectively . These models fit59 a different genetic variance for each environment and a 60 different genetic covariance for each pair of environments. 61 Model 2 comprises a single mean genetic value across environ- 62 ments for each genotype while the remaining models comprise 63 specific genetic values for each genotype by environment com- 64 bination. All models were fitted using ASReml-R (Butler et al. 65 2017), which obtains Residual Maximum Likelihood (REML) 66 estimates of the variance parameters and empirical Best Linear 67 Unbiased Predictions (BLUPs) of the random effects. Model fit 68 was assessed using the residual log-likelihood, AIC, and per- 69 centage of variance explained ( ve). The ve was calculated as 70 the variance explained by the correlated genetic values as a per- 71 centage of the total genetic variance following Tolhurst et al. 72 (2019). Models 2 and 4 were then used for CV and genome-wide 73 association, as described in the following. 74 We investigated how the different RMs impact the predictive 75 ability of the different linear mixed models. Due to generational 76 overlap in our data and a low number of individuals in the last 77 generations, we could not perform forward validation of the 78 genomic prediction of grain yield and have performed a five- 79 fold CV instead. To ensure a sufficiently structured and powered 80 CV , we used a subset of the phenotypic records corresponding 81 to the three years with the highest number of shared individuals 82 (2011-2013). The subset comprised 93 trials at Paso de la Laguna 83 with 327 indica and 330 japonica individuals in total. We used five- 84 fold CV , with the dataset partitioned into five sets containing 85 about 65 indica and 66 japonica individuals. Each set was used 86 as a validation set while the remaining four sets were used to 87 train the model for prediction. We tested five CV scenarios: 88 (i) CVI J→I J trained on and predicted both indica and japonica, 89 (ii) CVI→J trained on indica and predicted japonica, (iii) CVJ→I 90 trained on japonica and predicted indica, (iv) CVI→I trained on 91 and predicted indica, and (v) CVJ→J trained on and predicted 92 japonica. Each scenario was replicated 1,000 times. We computed 93 the predictive ability for each scenario as Pearson’s correlation 94 coefficient (r) between the predicted mean genetic values and 95 mean phenotypic values across environments. For model 4, r 96 was also computed between the predicted specific genetic values 97 and mean phenotypic values within environments. 98 We performed a genome-wide association study (GWAS) on 99 the complete data set, with standardised effects of the SNP sites 100 obtained following Gualdrón Duarte et al. (2014): 101 zi = ˜αip var( ˜αi) , (2) where ˜αi is the the BLUP of the ith SNP site effect and var( ˜αi) is 102 its variance. These components were obtained from the genetic 103 values parameterised by the SRM and BRM following Tolhurst 104 et al. (2019), as described in the Appendix. Note, however, since 105 the matrix D was not readily available for this study we were 106 limited to obtaining near-unique BLUPs and variances for the 107 BRM. The associated p-values from both RMs were displayed 108 using Manhattan plots, with significance assessed at p = 0.05 109 using Li and Ji (2005) correction that accounts for the effective 110 number of independent tests. For each model, we also calculated 111 Pearson’s correlation coefficient between the predicted SNP site 112 effects obtained with each RM. 113 Finally , we selected a region spanning two GWAS hits and examined how the effects of the SNP sites from model 2 with the BRM mapped to the branches of the tree sequence and cor- responding haplotype effects. A haplotype was defined as the .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 20, 2025. ; https://doi.org/10.1101/2025.01.14.633033doi: bioRxiv preprint Rebollo et al. 2025 5 combination of the allelic states of contiguous SNP sites. Haplo- type effects were estimated as a linear combination of their SNP site effects w⊤ j ˜α and standardised using: hj = w⊤ j ˜α q w⊤ j var( ˜α)wj , (3) where wj denotes a vector of 0s and 1s for allelic states for the jth1 haplotype and w⊤ j var( ˜α)wj the variance of the haplotype effect,2 with var( ˜α) being the variance-covariance matrix of ˜α. Note that3 our interpretation of the hj is limited due to the near-unique4 solutions of the ˜αi mentioned above.5 Results6 Overall, our results show that: (i) tree sequences captured the7 underlying biological signal from the genomic data; (ii) tree8 sequences encoded the genomic data more succinctly than the9 standard format; (iii) the site-based (SRM) and branch-based10 (BRM) relationship matrices were highly correlated, as were11 subsequent grain yield predictions; (iv) the BRM achieved the12 highest predictive ability when information from both indica and13 japonica was analyzed; and (iv) SNP site effects mapped onto14 branches of an ARG provide a way to estimate haplotype effects15 and give insight into the evolution of their value.16 The results are structured as follows. We initially demonstrate17 the usefulness of tree sequences to capture biological signal sum-18 marised with genealogical nearest neighbors (GNN) and local19 trees for regions of interest. We then summarise the ancestral20 alleles and tree sequences inference. The latter includes the21 summary of memory usage and compression, as well as the22 estimated age of ancestors, mutations, and SNP sites. Finally ,23 we present the results of quantitative genetic analyses based24 on tree sequences, including genomic prediction and GWAS.25 Here, we present the results for Ne = 1,500 when dating the tree26 sequences. The results for Ne = 23 and 150,000 are presented in27 the Supplementary Material.28 Our results show that tree sequences captured the underlying29 biological signal in the genomic data. The GNN captured the in-30 herent population structure and clearly separated the indica and31 japonica rice subspecies (Figure 1). Within subspecies, the GNN32 revealed a more pronounced structure for indica compared to33 japonica, possibly indicating a stronger family structure. Figure 134 also demonstrated admixed individuals between subspecies or35 potentially mislabeled individuals. We found a differential local36 structure by inspecting two interesting regions of the genome37 (Figure 2). Both regions had two SNP sites and two local trees38 (one SNP per tree), but the adjacent local trees were practically39 indistinguishable. The trees corresponding to the first SNP for40 both regions are shown in Figure 2 while the trees of the second41 SNP are shown in Supplementary Figure 2. The tree spanning42 the DST gene showed a clear and deep separation between in-43 dica and japonica (Figure 2A). Here, 98% of japonica individuals44 had the ancestral haplotype at both SNP sites, while 77% of45 indica individuals had a mutation at the first SNP and 19% at46 the second SNP (Supplementary Table 3). The tree covering the47 locus associated with NSB and NSSB showed that this region48 was segregating in both indica and japonica (Figure 2B). Here,49 27% of indica individuals and 80% of japonica individuals had50 the ancestral haplotype at both SNP sites while 71% of indica51 individuals had a mutation at the second SNP and 19% of japon-52 ica individuals had a mutation at the first SNP (Supplementary 53 Table 3). 54 We inferred ancestral alleles for about 80% of the loci in the 55 three genomic datasets (ALL, IND, and JAP; Table 1). Ancestral 56 alleles matched the major allele at 68% loci in the ALL dataset, 57 59% loci in the JAP dataset, and 52% loci in the IND dataset. 58 In less than 1% of all cases, the ancestral allele was neither the 59 major nor the minor allele. Some alleles were more common in 60 indica while others were more common in japonica. The distri- 61 bution of the alternative (non-ancestral) allele frequency within 62 each subspecies was U-shaped, with an excess of the lower alter- 63 native allele frequency for japonica (Figure 3A). However, when 64 combined, the distribution indicated differences in alternative 65 allele frequencies between the two subspecies (Figure 3B). Some 66 SNP sites had a low frequency in one subspecies and a high fre- 67 quency in the other, and vice-versa (Figure 3C). Similar results 68 were demonstrated by the distribution of minor allele frequen- 69 cies (Supplementary Figure 3). 70 The inferred tree sequences encoded the genomic data more 71 succinctly than the standard format. The combined tree se- 72 quence across chromosomes included 95,612 nodes, 693,524 73 edges, 31,925 local trees, 61,260 SNP sites, and 1,031,672 mu- 74 tations in total (Table 2). These statistics varied for individual 75 chromosomes, and were generally well-correlated with the chro- 76 mosome length. The combined file size of the dated trees across 77 chromosomes was 62 Mb while the standard VCF file was 228 78 Mb, i.e., almost four times larger. 79 The age distributions for nodes (ancestors), mutations, and 80 SNP sites (i.e., first mutation at each site) were heavily right- 81 skewed toward the present (Figure 4). For the nodes, this is 82 expected due to the coalescence of nodes backward in time. 83 Furthermore, since the dataset involved a breeding population, 84 we expect to observe rapid coalescence due to selection and 85 a small number of parents. Assuming Ne = 1,500, 70.4% of 86 nodes were estimated to be younger than 1,000 generations, 87 61.5% younger than 500 generations, and 41.2% younger than 88 100 generations. Conversely , 54.8% of mutations were estimated 89 to be younger than 2,500 generations, 30.4% were estimated 90 between 2,500 and 5,000 generations, and 14.8% were estimated 91 above 5,000 generations. All chromosomes showed a similar age 92 pattern for nodes, SNP sites, and mutations. The Ne used when 93 dating the tree sequences strongly affected the estimated age of 94 nodes, SNP sites, and mutations (Supplementary Figure 4). 95 Heatmaps of the SRM and BRM revealed a very similar pop- 96 ulation structure with elements of the matrices being highly 97 correlated, although on a different scale (Figure 5). The correla- 98 tion between the diagonal elements of the SRM and BRM was 99 0.98 across both subspecies (0.75 for indica and 0.89 for japonica) 100 while the correlation between the off-diagonal elements was 0.99 101 (1.00 for indica and 0.98 for japonica. The Ne used when dating 102 the tree sequences impacted the resulting BRM; with Ne = 23 103 producing significantly lower values (Supplementary Figure 5). 104 The predicted mean genetic values were highly correlated 105 between the SRM and BRM when the subspecies in the train- 106 ing and prediction sets were the same ( r = 0.96 − 0.99; Table 3). 107 However, the correlations were lower when the subspecies in the 108 training and prediction sets were different ( r = 0.74 − 0.88). The 109 equivalent correlations involving the PRM were much more var- 110 ied, e.g., the correlation was low ( CVI→J) and negative ( CVJ→I) 111 when the training and prediction sets were different. Finally , 112 note that the correlations involving the PRM were higher with 113 the BRM than with the SRM. 114 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 20, 2025. ; https://doi.org/10.1101/2025.01.14.633033doi: bioRxiv preprint 6 Genomic Prediction with Ancestral Recombination Graphs in rice Training in both subspecies together gave similar predic-1 tive abilities to training within a single subspecies only (Ta-2 ble 4). Models with the BRM had the highest predictive abil-3 ity for CVI J→I J and for the within subspecies scenarios CVI→I4 and CVJ→J. Overall, the predictive abilities were higher for in-5 dica than japonica. The highest predictive abilities ( ≥ 0.4) were6 achieved in CVI J→I J, CVI→I, and CVJ→J for all models and RMs.7 Conversely , the lowest predictive abilities were observed for the8 challenging scenario CVI→J, but not for CVJ→I where high pre-9 dictive abilities were achieved with the SRM and BRM, although10 negative predictive abilities were observed with the PRM. For11 training in one subspecies and predicting into the other ( CVI→J12 and CVJ→I), there was no consistently best RM. For training13 and predicting within the same subspecies ( CVI→I and CVJ→J),14 the BRM gave the highest predictive abilities. Overall, the PRM15 performed unexpectedly well except for CVJ→I. The BRM built16 with Ne of 23, 1500, and 150,000 were compared for a favorable17 scenario (CVI J→I J) and an unfavorable scenario ( CVI→J ), and18 predictive abilities with Ne = 1500 and 150,000 were the highest19 and very similar (Supplementary Table 7).20 The predicted SNP site effects from models 2 and 4 were21 highly correlated between the SRM and BRM (Table 5). The22 number of significant SNP site effects for each model and RM23 was 13 for model 2 and 14 for model 4 with the SRM and 1024 for model 2 and 5 for model 4 with the BRM (Figure 6). The25 standardised SNP site effects for these models are presented26 in Supplementary Figure 6. A total of 27 SNP site effects were27 significant across models and RMs. No SNP site effect was28 significant across all model and RM combinations. Specifically ,29 one SNP site effect was significant across 3 combinations, 1330 SNP site effects were significant across two combinations, and31 the remaining 13 were unique to a single combination.32 We used the zoom-in region to study how the estimated SNP33 site effects map onto ARG branches and generate haplotype34 effects. We studied the region between GWAS hits at positions35 29,476,724 and 29,561,694 of chromosome 6, which spans 84,97036 bp, corresponding to ∼ 0.3cM. We observed 15 SNP sites in 837 local trees within this region (Supplementary Table 4). Of the 838 trees, 2 contained 4 SNP sites each, 1 tree contained 2 SNP sites,39 and the remaining 5 contained 1 SNP each. The frequency of40 SNP mutations was generally higher in indica than in japonica.41 The standardised SNP site effect was between -3.55 and 3.5242 (Supplementary Table 4). There were a total of 61 haplotypes in43 the region. The haplotype effect of each local tree ranged from44 -3.55 to 5.77 for the trees with 4 SNP sites, from 0.00 to 2.58 for45 the tree with 2 SNP sites, and from -2.78 to 3.28 for the trees with46 1 SNP (Supplementary Table 5). Finally , the absolute value of47 the haplotype effect of the 10 most frequent haplotypes ranged48 from 0.00 to 3.93 (Supplementary Table 6).49 Discussion50 This study adds to the growing body of literature on quantitative51 genetic analysis across populations (Martin et al. 2017; Ding et al.52 2023; Mester et al. 2023; Warburton et al. 2023; Zhang et al. 2023b)53 and the applications and benefits of ARGs (Kelleher et al. 2019;54 Wohns et al. 2022; Fan et al. 2022; Brandt et al. 2022; Link et al.55 2023; Zhang et al. 2023a; Brandt et al. 2024; Lewanski et al. 2024;56 Sun et al. 2024; Gunnarsson et al. 2024; Zhu et al. 2024; Nielsen57 et al. 2024). We evaluated the use of an ARG to capture the ge-58 netic history and improve quantitative genetic modelling in a59 rice breeding dataset with indica and japonica subspecies. In the60 following, we discuss the usefulness of an ARG encoded as a61 tree sequence for capturing biological features of the analyzed 62 population and for efficiently storing genomic data. We also 63 discuss how tree sequences can capture relationships between 64 the genomes of individuals, making them highly relevant to ge- 65 nomic prediction and GWAS. Finally , we discuss the limitations 66 of the study , particularly the challenges with inferring ARGs and 67 the impact of estimating Ne on dating. Overall, our results high- 68 light the potential of ARGs as an emerging tool for quantitative 69 genetic analyses and their application in selective breeding. 70 T ree sequences for capturing biological signal and data 71 compression 72 We corroborated that tree sequences are useful in portraying 73 population structure and are efficient in compressing genomic 74 data. Our GNN analysis successfully captured the underly- 75 ing population structure and corresponding biological signals; 76 clearly separating the indica and japonica rice subspecies. The 77 clear within-subspecies structure can be explained by the rice 78 breeding program, which uses a few key individuals as par- 79 ents for many crosses. The GNN analysis also proved useful 80 for identifying mislabeled individuals and determining the sub- 81 species of unlabelled individuals, increasing their usefulness in 82 future genetic analyses. We found differential local structure by 83 inspecting individual trees in two genome regions of interest 84 concordant with previous studies (Bai et al. 2016). Although 85 there were only two SNP sites in each region, we observed a 86 deep split at the DST locus between indica and japonica with a 87 low frequency of both mutations in indica. Conversely , we ob- 88 served that the locus associated with NSB and NSSB segregated 89 in both subspecies. Although there is insufficient data to con- 90 clude whether the genomic differentiation observed in the trees 91 is associated with any phenotypic differentiation, the analyses 92 performed here provide greater understanding of the origin and 93 evolution of known QTLs. Such genealogical insight can also 94 aid in differentiating selection and causality from demography 95 responsible for spurious GWAS hits. 96 The tree sequence format has been shown to be extremely 97 beneficial for dense and whole-genome sequence data (Kelleher 98 et al. 2019), and here it compressed the data almost 4 times more 99 than a standard VCF file. Despite this, further work is needed 100 to improve ARG inference in the rice dataset. For example, we 101 inferred 1,031,672 mutations across 61,260 SNP sites, producing 102 an average of ∼ 17 mutations per site, which is an unexpected 103 result. We followed the approach in Wohns et al. (2022), who 104 observed 5,773,816 mutations across 2,090,401 variable SNP sites, 105 producing an average of ∼ 3 mutations per site. The key differ- 106 ence is that they used whole-genome sequencing data, whereas 107 we used reduced-representation GBS data. The high number of 108 mutations per site here may be attributed to the chosen ARG in- 109 ference parameters; sequencing, imputation and phasing errors; 110 insufficient number of SNP sites; errors in the inferred ancestral 111 allele; or the ARG inference method. While some imputation 112 and phasing errors are expected, the rice individuals are highly 113 inbred so this an unlikely issue, even if the default parameters 114 in Beagle are not optimal for inbred crops such as rice (Niehoff 115 et al. 2022). 116 An alternative inference method may overcome the issues 117 above, e.g., the infinite-sites model, which allows only one mu- 118 tation per site and thus produces a higher number of inferred 119 recombination events and local trees (Wohns et al. 2022; Kelleher 120 et al. 2019). Further research is also warranted to optimise the 121 designation of mutation and recombination events (Wohns et al. 122 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 20, 2025. ; https://doi.org/10.1101/2025.01.14.633033doi: bioRxiv preprint Rebollo et al. 2025 7 2022). Such research could use ARGneedle (Zhang et al. 2023a)1 instead of tsinfer (Kelleher et al. 2019, 2021b), to also address2 SNP site ascertainment bias.3 T ree sequences for quantitative genetic analyses4 We showed that tree sequences are useful for capturing relation-5 ships between the genomes of individuals and can be utilised6 for genomic prediction and GWAS. The SRM and BRM revealed7 very similar underlying population structure and were highly8 correlated, as expected given the duality between these two9 statistics (Ralph 2019; Ralph et al. 2020). Consequently , we ob-10 served a high correlation between the predicted mean genetic11 values obtained with the SRM and BRM for all CV scenarios,12 particularly when the subspecies in the training and prediction13 sets were the same. The equivalent correlations involving the14 PRM were much lower, and even negative when the subspecies15 in the training and prediction sets were different. This is because16 predicting from indica to japonica and vice-versa is very difficult17 with the PRM due to the lack of deep pedigree and the large18 number of meiosis events between the subspecies. It is worth19 noting, however, that the correlations involving the PRM were20 higher with the BRM than with the SRM for all CV scenarios.21 This may be attributed to the fact that the BRM and PRM both22 capture identity by descent information; albeit the BRM captures23 realised identity by descent from the roots of an ARG and the24 PRM captures expected identity by descent from the pedigree25 founders.26 All models achieved a high predictive ability when the sub-27 species in the training and prediction sets were the same. This is28 expected since the accuracy of genomic prediction is a function29 of genetic distance between training and prediction individuals30 (Hickey et al. 2014; Lorenz and Smith 2015; Scutari et al. 2016;31 Ding et al. 2023). In general, we observed that the BRM achieved32 the highest predictive ability for these scenarios, highlighting33 the benefit of ARGs in jointly modelling identity by descent and34 identity by state relationships. Modest improvements were ob-35 served over the SRM, likely due to the limited dataset size, much36 larger number of parameters estimated (693,524 ARG edges com-37 pared to 61,260 SNP sites), and potential errors arising in the38 underlying ARG.39 Training in both subspecies together gave similar predictive40 abilities to training within a single subspecies only . This is con-41 trary to a previous study involving a subset of our data Berro42 et al. (2019), and could be explained by the larger number of43 SNP sites used here (61,260 SNP sites compared to 15,545). We44 observed that indica had higher predictive ability than japonica45 for these scenarios, which is expected due to large groups of46 indica relatives within the breeding program (Hickey et al. 2014;47 Scutari et al. 2016; Ding et al. 2023). Finally , training in one sub-48 species and predicting into the other subspecies produced the49 lowest predictive abilities of all scenarios. This is in agreement50 with previous empirical studies that reported limited predic-51 tive ability between populations, let alone between different52 subspecies like indica and japonica rice. The results may there-53 fore suggest the presence of a different genetic architecture for54 complex traits between the subspecies (Zhao et al. 2011), dif-55 ferent linkage-disequilibrium between typed loci and QTL, or56 additional (non-additive) genetic effects.57 We observed correspondence between the predicted SNP site58 effects obtained with the SRM and BRM, indicating consistency59 across models and RMs. Specifically , the significant SNP site60 effects were spread across different model and RM combina-61 tions. Approximately half of the significant SNP site effects were 62 unique to specific model and RM combinations, while the other 63 half were shared across two combinations. Focusing on a spe- 64 cific region of chromosome 6 between two GWAS hits provided 65 additional insights into haplotype effects, revealing variation in 66 SNP and haplotype effects across trees. Such insights will be en- 67 hanced in future studies by the presence of known relationship 68 information between SNP sites, allowing unique predictions of 69 site and haplotype effects to be obtained. 70 Considerations 71 Plant breeding datasets are challenged by intense selection from 72 one cycle to the next. We hypothesize that a forward validation 73 scheme could highlight the benefit of more precisely capturing 74 population and family structure with ARGs compared to the 75 five-fold CV used here. Furthermore, the size and population 76 structure of our dataset may be the limiting factor as to why the 77 genomic RMs did not produce much higher predictive abilities 78 than the PRM for the more favourable scenarios CVI J→I J, CVI→I, 79 and CVJ→J (Daetwyler et al. 2010; Lorenz and Smith 2015; Papin 80 et al. 2024). However, a large portion of the individuals in the 81 dataset did not have genomic data, limiting the study to only 82 those with genomic data available. For these cases, the integra- 83 tion of pedigree and genomic data has been proven successful 84 (Aguilar et al. 2010), but future work is needed to integrate them 85 in a single-step version using the BRM. 86 Depending on the method and software implementation, 87 ARG inference and dating require a high density of SNP sites, 88 knowledge about Ne, ancestral alleles, and mutation and re- 89 combination rate. Although rice’s mutation and recombination 90 rate estimates are available (Yang et al. 2015; Si et al. 2015), the 91 ancestral alleles and Ne for our population were not. For 25% of 92 the SNP sites, we were not able to infer the ancestral allele and 93 assumed the major allele as the ancestral. This likely introduced 94 errors, since the chosen ancestral allele was likely to be driven 95 by the higher number of japonica individuals in the combined 96 dataset. 97 The choice of Ne is not trivial in agricultural species, since 98 it changes significantly over time due to domestication, bottle- 99 necks, and selective breeding. For rice, we found very contrast- 100 ing estimates related to different time periods; current 23, recent 101 1500, and ancient 150,000. The Ne assumption strongly impacted 102 the dating of trees, but not the predictive ability of the different 103 models. Finally , we did not use high-density whole-genome 104 sequence data but instead used GBS data with 61,260 SNP sites 105 which can limit the ARG inference. Note, however, it has been 106 shown that accurate ARG inference can be achieved for subsets 107 of whole-genome sequence data, such as SNP arrays, when Ne 108 is small (Fan et al. 2022). 109 Future directions 110 There are many new opportunities to expand this research. 111 Firstly , more research is warranted to improve the inference 112 and dating of ARGs, and to understand the inference behaviour 113 for multiple species. Secondly , large-scale datasets are required 114 to achieve sufficient power to estimate the very large number of 115 edge effects, and efficient algorithms are currently underway to 116 compute the BRM for such datasets (P . Ralph, personal commu- 117 nication, August 30th, 2024). Finally , ARGs have the potential to 118 improve genomic prediction and GWAS across diverse popula- 119 tions as they can capture mutations at untyped SNP sites on the 120 inferred edges between genomes (Selle et al. 2021; Harris 2023; 121 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 20, 2025. ; https://doi.org/10.1101/2025.01.14.633033doi: bioRxiv preprint 8 Genomic Prediction with Ancestral Recombination Graphs in rice Link et al. 2023; Zhang et al. 2023a; Brandt et al. 2024; Lewanski1 et al. 2024; Gunnarsson et al. 2024; Zhu et al. 2024), and may2 also capture novel epistatic variations responsible for changing3 mutation effects over time (Mathieson 2021; Park et al. 2022).4 Conclusions5 We have shown that: (i) an ARG encoded as a tree sequence6 effectively captured population structure and underlying bio-7 logical signal in indica and japonica subspecies in the Uruguayan8 National Rice Breeding Program. (ii) Tree sequences can be used9 for genomic prediction and GWAS, and have great potential in10 larger datasets by estimating mutation effects in different genetic11 backgrounds. (iii) In terms of GWAS, inspecting the ancestry of12 haplotypes and their values provides insight into the history of13 mutations and resulting haplotype differences. (iv) There are14 still challenges with accurately inferring ARGs. In conclusion,15 our results highlight the potential of ARGs as an emerging tool16 for the quantitative genetic analysis of diverse populations.17 Data availability18 The dataset, including phenotypic and genomic data, is publicly19 available at the Dryad repository (Rebollo et al. 2023b) with20 anonymized line identification. Analysis scripts and data are21 available at https://github.com/HighlanderLab/irebollo_rice_tree.22 Contributions23 IR curated and analyzed the data, prepared the results, drafted,24 and edited the manuscript. DT guided the statistical analy-25 sis, derived and wrote the appendix for backsolving SNP site26 effects from genetic values, contributed to the interpretation27 of the results, and edited the manuscript. JO contributed to28 tree sequences inference, interpretation of results, and edited29 the manuscript. JR led the engagement with the rice breeding30 program, contributed to data generation and curation, inter-31 pretation of results, and edited the manuscript. GG led the32 study , contributed to the interpretation of results, and edited the33 manuscript.34 Acknowledgments35 The authors acknowledge INIA’s rice breeding team: Pedro36 Blanco, Fernando Pérez de Vida, Federico Molina, and former37 and current field and laboratory staff. We also acknowledge38 Brieuc Lehmann, Jerome Kelleher, and Peter Ralph for tskit39 implementation of the branch-based relationship matrix (BRM),40 and Yan Wong for discussions about tree sequences.41 Funding42 The authors acknowledge funding from Instituto Nacional de43 Investigación Agropecuaria, Uruguay (Projects AZ35, AZ13,44 and fellowship to IR), Agencia Nacional de Investigación e In-45 novación, Uruguay (Project FSDA_1_2018_1_154120), Comité46 Académico de Posgrado (fellowship to IR). The authors also47 acknowledge funding from the Edinburgh Innovations Fellow-48 ship to DT, the support from the Slovenian Research Agency’s49 research program P4-0133 to JO, and the BBSRC ISP grant to50 The Roslin Institute (BBS/E/D/30002275, BBS/E/RL/230001A,51 BBS/E/RL/230001C), BBSRC projects BB/R019940/1 and52 BB/T014067/1, and The University of Edinburgh.53 Appendix 54 This appendix derives an approach to backsolve for SNP site 55 effects from genetic values parameterised by a site-based (SRM) 56 or branch-based (BRM) relationship matrix. Assume the vector 57 of random genetic values in Eq. 1 are modelled as u = Mα, 58 where M is a centered SNP genotype matrix and α is a vector 59 of SNP site effects, often referred to as allele substitution effects 60 when parameterised by the SRM. 61 It is assumed that: 62   u α   ∼ N     0 0   , σ2 α   G MD DM⊤ D     , where σ2α is the variance parameter of the SNP site effects, 63 G = MDM⊤ is the relationship matrix between individuals 64 and D is the relationship matrix between SNP site effects. Typi- 65 cally , the matrixD is assumed to be diagonal (VanRaden 2008), 66 like for the SRM where D = I/c and c is the sum of the expected 67 heterozygosities across SNP sites. In terms of the BRM, however, 68 D was unknown and assumed to be non-diagonal, with relation- 69 ships between sites induced through the tree sequence structure. 70 In this case, we approximated D; as: 71 D ≈ M−GM−⊤, (4) where M− is the Moore-Penrose generalised inverse of M, ob- 72 tained such that D satisfies G ≡ MDM⊤. We have found that D 73 obtained in this manner represents a near-unique solution, and 74 have made brief interpretation in text to reflect this limitation. 75 Let ˜u denote the vector of BLUPs of the genetic values param- 76 eterised by the site- or branch-based RM. Following Tolhurst 77 et al. (2019), the BLUPs of the SNP site effects can be obtained 78 from ˜u as: 79 ˜α = DM⊤G−1 ˜u, (5) with corresponding variance matrix obtained as: 80 var( ˜α) = σ2 α DM⊤G−1MD − DM⊤G−1PEV( ˜u)G−1MD, (6) where PEV( ˜u) is the prediction error variance matrix of the 81 predicted genetic values, routinely supplied by mixed model 82 software. The test statistics in Eq. 2 are then constructed using 83 the components in Eqs. 5 and 6. 84 Conflicts of interest 85 The authors declare no conflict of interest. 86 Literature cited 87 Aguilar I, Misztal I, Johnson D, Legarra A, Tsuruta S, Lawlor T. 88 2010. Hot topic: A unified approach to utilize phenotypic, full 89 pedigree, and genomic information for genetic evaluation of 90 Holstein final score1. Journal of Dairy Science. 93:743–752. 91 Aguilar I, Misztal I, Tsuruta S, Legarra A. 2014. PREGSF90 – 92 POSTGSF90 : Computational Tools for the Implementation 93 of Single-step Genomic Selection and Genome-wide ... Pro- 94 ceedings of the 10th World Congress of Genetics Applied to 95 Livestock Production. V olume Met:1–3. 96 Bai X, Zhao H, Huang Y , Xie W, Han Z, Zhang B, Guo Z, Yang L, 97 Dong H, Xue W et al. 2016. Genome-Wide Association Analy- 98 sis Reveals Different Genetic Control in Panicle Architecture 99 Between Indica and Japonica Rice . The Plant Genome. 9:1–10. 100 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 20, 2025. ; https://doi.org/10.1101/2025.01.14.633033doi: bioRxiv preprint Rebollo et al. 2025 9 Batello C, Queheille N, Uraga R, Gonnet D, Hill M, Roel A,1 Cantou G, Martínez I, Pippolo D. 2013. Guia de buenas prácticas2 en el cultivo de arroz en Uruguay.. Asociación de Cultivadores3 de Arroz, ACA. Uruguay .4 Berro I, Lado B, Nalin RS, Quincke M, Gutiérrez L. 2019. Training5 Population Optimization for Genomic Selection. The Plant6 Genome. 12:1–14.7 Bradbury PJ, Zhang Z, Kroon DE, Casstevens TM, Ramdoss Y ,8 Buckler ES. 2007. TASSEL: Software for association mapping9 of complex traits in diverse samples. Bioinformatics. 23:2633–10 2635.11 Brandt DYC, Huber CD, Chiang CWK, Ortega-Del V ecchyo D.12 2024. The Promise of Inferring the Past using the Ancestral13 Recombination Graph (ARG). Genome Biology and Evolution.14 p. evae005.15 Brandt DYC, Wei X, Deng Y , Vaughn AH, Nielsen R. 2022. Eval-16 uation of methods for estimating coalescence times using an-17 cestral recombination graphs. Genetics. 221.18 Browning BL, Zhou Y , Browning SR. 2018. A One-Penny Im-19 puted Genome from Next-Generation Reference Panels. The20 American Journal of Human Genetics. 103:338–348.21 Butler DG, Cullis BR, Gilmour AR, Gogel BJ, Thompson R. 2017.22 ASReml-R Reference Manual V ersion 4. ASReml-R Reference23 Manual. p. 176.24 Caicedo AL, Williamson SH, Hernandez RD, Boyko A, Fledel-25 Alon A, York TL, Polato NR, Olsen KM, Nielsen R, McCouch26 SR et al. 2007. Genome-wide patterns of nucleotide polymor-27 phism in domesticated rice. PLoS Genetics. 3:1745–1756.28 Corcoran P , Barton HJ. 2021. Wgabed. https://henryjuho.github.io/29 WGAbed/.30 Daetwyler HD, Hickey JM, Henshall JM, Dominik S, Gredler31 B, Van Der Werf JH, Hayes BJ. 2010. Accuracy of estimated32 genomic breeding values for wool and meat traits in a multi-33 breed sheep population. Animal Production Science. 50:1004–34 1010.35 de los Campos G, Hickey JM, Pong-Wong R, Daetwyler HD,36 Calus MP. 2013. Whole-genome regression and prediction37

applied to plant and animal breeding. Genetics.38 193:327–345.39 Ding Y , Hou K, Xu Z, Pimplaskar A, Petter E, Boulier K, Privé40 F, Vilhjálmsson BJ, Olde Loohuis LM, Pasaniuc B. 2023. Poly-41 genic scoring accuracy varies across the genetic ancestry con-42 tinuum. Nature. 618:774–781.43 Elshire RJ, Glaubitz JC, Sun Q, Poland JA, Kawamoto K, Buck-44 ler ES, Mitchell SE. 2011. A robust, simple genotyping-by-45 sequencing (GBS) approach for high diversity species. PloS46 one. 6:e19379.47 Endelman JB. 2011. Ridge Regression and Other Kernels for48 Genomic Selection with R Package rrBLUP. The Plant Genome.49 4:250–255.50 Fan C, Mancuso N, Chiang CW. 2022. A genealogical estimate51 of genetic relationships. The American Journal of Human52 Genetics. 109:812–824.53 Garris AJ, Tai TH, Coburn J, Kresovich S, McCouch S. 2005.54 Genetic structure and diversity in Oryza sativa L. Genetics.55 169:1631–1638.56 Glaszmann JC. 1987. Isozymes and classification of Asian rice57 varieties. Theoretical and Applied Genetics. 74:21–30.58 Griffiths RC, Marjoram P. 1996. Ancestral inference from sam-59 ples of dna sequences with recombination. Journal of compu-60 tational biology . 3:479–502.61 Gualdrón Duarte JL, Cantet RJ, Bates RO, Ernst CW, Raney62 NE, Steibel JP. 2014. Rapid screening for phenotype-genotype 63 associations by linear transformations of genomic evaluations. 64 BMC Bioinformatics. 15:246. 65 Gunnarsson AF, Zhu J, Zhang BC, Tsangalidou Z, Allmont A, 66 Palamara PF. 2024. A scalable approach for genome-wide 67 inference of ancestral recombination graphs. bioRxiv . . 68 Haller BC, Galloway J, Kelleher J, Messer PW, Ralph PL. 2019. 69 Tree-sequence recording in SLiM opens new horizons for 70 forward-time simulation of whole genomes. Molecular Ecol- 71 ogy Resources. 19:552–566. 72 Harris K. 2023. Using enormous genealogies to map causal vari- 73 ants in space and time. Nature Genetics. 55:730–731. 74 Henderson CR. 1984. Applications of linear models in animal breed- 75 ing. University of Guelph. Guelph. 76 Hickey JM, Dreisigacker S, Crossa J, Hearne S, Babu R, Prasanna 77 BM, Grondona M, Zambelli A, Windhausen VS, Mathews K 78 et al. 2014. Evaluation of Genomic Selection Training Popu- 79 lation Designs and Genotyping Strategies in Plant Breeding 80 Programs Using Simulation. Crop Science. 54:1476–1488. 81 Huang X, Kurata N, Wei X, Wang ZX, Wang A, Zhao Q, Zhao Y , 82 Liu K, Lu H, Li W et al. 2012. A map of rice genome variation 83 reveals the origin of cultivated rice. Nature. 490:497–501. 84 Kawahara Y , de la Bastide M, Hamilton JP , Kanamori H, Mccom- 85 bie WR, Ouyang S, Schwartz DC, Tanaka T, Wu J, Zhou S et al. 86 2013. Improvement of the oryza sativa nipponbare reference 87 genome using next generation sequence and optical map data. 88 Rice. 6:3–10. 89 Keightley PD, Jackson BC. 2018. Inferring the probability of the 90 derived vs. The ancestral allelic state at a polymorphic site. 91 Genetics. 209:897–906. 92 Kelleher J, Etheridge AM, McV ean G. 2016. Efficient coalescent 93 simulation and genealogical analysis for large sample sizes. 94 PLoS computational biology . 12:e1004842. 95 Kelleher J, Jeffery B, Wong Y , Ralph P , Tsambos G, Goldstein D, 96 Lehmann B, Thornton KR, Wohns AW, pyup.io bot et al. 2021a. 97 tskit-dev/tskit: Python 0.3.7. 98 Kelleher J, Thornton KR, Ashander J, Ralph PL. 2018. Efficient 99 pedigree recording for fast population genetics simulation. 100 PLoS Computational Biology . 14:1–21. 101 Kelleher J, Wong Y , Jeffery B, Wohns AW, Rebollo I, F DR. 2021b. 102 tskit-dev/tsinfer: 0.2.1. 103 Kelleher J, Wong Y , Wohns AW, Fadil C, Albers PK, McV ean G. 104 2019. Inferring whole-genome histories in large population 105 datasets. Nature Genetics. 51:1330–1338. 106 Lewanski AL, Grundler MC, Bradburd GS. 2024. The era of the 107 arg: An introduction to ancestral recombination graphs and 108 their significance in empirical evolutionary genomics. PLOS 109 Genetics. 20:1–24. 110 Li H, Durbin R. 2009. Fast and accurate short read alignment 111 with Burrows–Wheeler transform. Bioinformatics. 25:1754– 112 1760. 113 Li J, Ji L. 2005. Adjusting multiple testing in multilocus analyses 114 using the eigenvalues of a correlation matrix. Heredity . 95:221– 115 227. 116 Link V , Schraiber JG, Fan C, Dinh B, Mancuso N, Chiang CWK, 117 Edge MD. 2023. Tree-based QTL mapping with expected local 118 genetic relatedness matrices. The American Journal of Human 119 Genetics. 110:2077–2091. 120 Liu D, Yin H, Li T, Wang L, Lu S, Yu H, Sun X, Zhang T, Zhao Q, 121 ZHAO Y. 2024. The origin and evolution of cultivated rice and 122 genomic signatures of heterosis for yield traits in super-hybrid 123 rice. bioRxiv . . 124 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 20, 2025. ; https://doi.org/10.1101/2025.01.14.633033doi: bioRxiv preprint 10 Genomic Prediction with Ancestral Recombination Graphs in rice Lorenz AJ, Smith KP. 2015. Adding genetically distant indi-1 viduals to training populations reduces genomic prediction2 accuracy in Barley . Crop Science. 55:2657–2667.3 Martin AR, Gignoux CR, Walters RK, Wojcik GL, Neale BM,4 Gravel S, Daly MJ, Bustamante CD, Kenny EE. 2017. Human5 demographic history impacts genetic risk prediction across di-6 verse populations. The American Journal of Human Genetics.7 100:635–649.8 Martínez CP , Torres EA, Chatel M, Mosquera G, Duitama J,9 Ishitani M, Selvaraj M, Dedicova B, Tohme J, Grenier C et al.10 2014. Rice breeding in latin america, In: , John Wiley & Sons,11 Ltd. chapter 5, pp. 187–278.12 Mathieson I. 2021. The omnigenic model and polygenic pre-13 diction of complex traits. The American Journal of Human14 Genetics. 108:1558–1563.15 Mester R, Hou K, Ding Y , Meeks G, Burch KS, Bhattacharya A,16 Henn BM, Pasaniuc B. 2023. Impact of cross-ancestry genetic17 architecture on gwas in admixed populations. bioRxiv . .18 Mrode. 2014. Linear Models for the Prediction of Animal Breeding19 Values. CABI. Boston. third edition.20 Niehoff T, Pook T, Gholami M, Beissinger T. 2022. Imputation of21 low-density marker chip data in plant breeding: Evaluation of22

based on sugar beet. The Plant Genome. 15:e20257.23 Nielsen R, Vaughn AH, Deng Y. 2024. Inference and applications24 of ancestral recombination graphs. Nature Reviews Genetics. .25 Papin V , Gorjanc G, Pocrnic I, Bouffier L, Sanchez L. 2024. Un-26 locking genome-based prediction and selection in conifers:27 the key role of within-family prediction accuracy illustrated in28 maritime pine (Pinus pinaster ait.). Annals of Forest Science.29 81.30 Park Y , Metzger BPH, Thornton JW. 2022. Epistatic drift causes31 gradual decay of predictability in protein evolution. Science.32 376:823–830.33 Pérez-Enciso M. 1995. Use of the uncertain relationship matrix to34 compute effective population size. Journal of Animal Breeding35 and Genetics. 112:327–332.36 Ralph P , Thornton K, Kelleher J. 2020. Efficiently summarizing37 relationships in large samples: A general duality between38 statistics of genealogies and genomes. Genetics. 215:779–797.39 Ralph PL. 2019. An empirical approach to demographic in-40 ference with genomic data. Theoretical population biology .41 127:91–101.42 Rebollo I, Aguilar I, Pérez de Vida F, Molina F, Gutiérrez L,43 Rosas JE. 2023a. Genotype by environment interaction charac-44 terization and its modeling with random regression to climatic45 variables in two rice breeding populations. Crop Science. .46 Rebollo I, Rosas J, Aguilar I. 2020. Efficient com- puta-47 tion of the additive relationship matrix and its inverse in48 self-breeding individuals. https://icqg6.org/wp-content/uploads/49 2021/04/ines-rebollo-EfficientComputation.mp4.50 Rebollo I, Scheffel S, Blanco P , Molina F, Martínez S, Carracelas51 G, Pérez de Vida F, Rosas JE. 2023b. Instituto nacional de inves-52 tigación agropecuaria (inia) rice breeding program historical53 dataset. [Dataset].54 Rutkoski JE. 2019. Estimation of realized rates of genetic gain55 and indicators for breeding program assessment. Crop Science.56 59:981–993.57 Santiago E, Novo I, Pardiñas AF, Saura M, Wang J, Caballero58 A. 2020. Recent Demographic History Inferred by High-59 Resolution Analysis of Linkage Disequilibrium. Molecular60 Biology and Evolution. 37:3642–3653.61 Scutari M, Mackay I, Balding D. 2016. Using Genetic Distance62 to Infer the Accuracy of Genomic Prediction. PLOS Genetics. 63 12:e1006288. 64 Selle ML, Steinsland I, Lindgren F, Brajkovic V , Cubric-Curik V , 65 Gorjanc G. 2021. Hierarchical Modelling of Haplotype Effects 66 on a Phylogeny. Frontiers in Genetics. 11:1–16. 67 Si W, Yuan Y , Huang J, Zhang X, Zhang Y , Zhang Y , Tian D, Wang 68 C, Yang Y , Yang S. 2015. Widely distributed hot and cold spots 69 in meiotic recombination as shown by the sequencing of rice 70 f2 plants. New Phytologist. 206:1491–1502. 71 Speidel L, Forest M, Shi S, Myers SR. 2019. A method for genome- 72 wide genealogy estimation for thousands of samples. Nature 73 Genetics. 51:1321–1329. 74 Sun Q, Rowland BT, Chen J, Mikhaylova A V , Avery C, Peters 75 U, Lundin J, Matise T, Buyske S, Tao R et al. 2024. Improving 76 polygenic risk prediction in admixed populations by explic- 77 itly modeling ancestral-differential effects via gaudi. Nature 78 Communications. 15:1016. 79 Tolhurst DJ, Gaynor RC, Gardunia B, Hickey JM, Gorjanc G. 80 2022. Genomic selection using random regressions on known 81 and latent environmental covariates. Theoretical and Applied 82 Genetics. 135:3393–3415. 83 Tolhurst DJ, Mathews KL, Smith AB, Cullis BR. 2019. Genomic 84 selection in multi-environment plant breeding trials using a 85 factor analytic linear mixed model. Journal of Animal Breed- 86 ing and Genetics. 136:279–300. 87 VanRaden PM. 2008. Efficient Methods to Compute Genomic 88 Predictions. Journal of Dairy Science. 91:4414–4423. 89 Virtanen P , Gommers R, Oliphant TE, Haberland M, Reddy T, 90 Cournapeau D, Burovski E, Peterson P , Weckesser W, Bright J 91 et al. 2020. SciPy 1.0: Fundamental Algorithms for Scientific 92 Computing in Python. Nature Methods. 17:261–272. 93 Warburton CL, Costilla R, Engle BN, Moore SS, Corbet NJ, 94 Fordyce G, McGowan MR, Burns BM, Hayes BJ. 2023. Con- 95 currently mapping quantitative trait loci associations from 96 multiple subspecies within hybrid populations. Heredity . pp. 97 1–11. 98 Wohns AW, Wong Y , Jeffery B, Akbari A, Mallick S, Pinhasi R, 99 Patterson N, Reich D, Kelleher J, McV ean G. 2022. A unified 100 genealogy of modern and ancient genomes. Science. 375. 101 Wong Y , Ignatieva A, Koskela J, Gorjanc G, Wohns AW, Kelleher 102 J. 2024. A general and efficient representation of ancestral 103 recombination graphs. Genetics. 228:iyae100. 104 Yang S, Wang L, Huang J, Zhang X, Yuan Y , Chen JQ, Hurst 105 LD, Tian D. 2015. Parent-progeny sequencing indicates higher 106 mutation rates in heterozygotes. Nature. 523:463–467. 107 Yates A, Allen J, Amode RM, Azov AG, Barba M, Becerra A, 108 Bhai J, Campbell L, Carbajo Martinez M, Chakiachvili M et al. 109 2021. Ensembl Genomes 2022: an expanding genome resource 110 for non-vertebrates. Nucleic Acids Research. 50:D996–D1003. 111 Young AI. 2022. Discovering missing heritability in whole- 112 genome sequencing data. Nature Genetics. 54:224–226. 113 Yu J, Pressoir G, Briggs WH, Bi IV , Yamasaki M, Doebley JF, 114 McMullen MD, Gaut BS, Nielsen DM, Holland JB et al. 2006. 115 A unified mixed-model method for association mapping that 116 accounts for multiple levels of relatedness. Nature Genetics. 117 38:203–208. 118 Zhang BC, Biddanda A, Gunnarsson ÁF, Cooper F, Palamara 119 PF. 2023a. Biobank-scale inference of ancestral recombination 120 graphs enables genealogical analysis of complex traits. Nature 121 Genetics. 55:768–776. 122 Zhang H, Zhan J, Jin J, Zhang J, Lu W, Zhao R, Ahearn TU, 123 Yu Z, O’Connell J, Jiang Y et al. 2023b. A new method for 124 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 20, 2025. ; https://doi.org/10.1101/2025.01.14.633033doi: bioRxiv preprint Rebollo et al. 2025 11 multiancestry polygenic prediction improves performance1 across diverse populations. Nature Genetics. pp. 1–12.2 Zhao K, Tung CW, Eizenga GC, Wright MH, Ali ML, Price AH,3 Norton GJ, Islam MR, Reynolds A, Mezey J et al. 2011. Genome-4 wide association mapping reveals a rich genetic architecture5 of complex traits in Oryza sativa. Nature Communications.6 2:467.7 Zhu J, Kalantzis G, Pazokitoroudi A, Gunnarsson ÁF, Loya H,8 Chen H, Sankararaman S, Palamara PF. 2024. Fast variance9 component analysis using large-scale ancestral recombination10 graphs. bioRxiv . .11 Zhu Q, Ge S. 2005. Phylogenetic relationships among A-genome12 species of the genus Oryza revealed by intron sequences of13 four nuclear genes. New Phytologist. 167:249–265.14 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 20, 2025. ; https://doi.org/10.1101/2025.01.14.633033doi: bioRxiv preprint 12 Genomic Prediction with Ancestral Recombination Graphs in rice Table 1 Characterization of the ancestral allele inference for the genomic datasets with all individuals from both subspecies (ALL) and from each subspecies separately (IND for indica and JAP for japonica). Percentages of the total number of ancestral alleles in- ferred are presented in parentheses. Genomic dataset T otal sites Alleles inferred Ancestral = Major Ancestral = Minor Ancestral ̸= Minor ̸= Major ALL 61,260 49,518 (80.8%) 33,485 (67.6%) 15,802 (31.9%) 231 (0.5%) IND 50,854 40,891 (80.4%) 21,253 (52.0%) 19,363 (47.3%) 275 (0.7%) JAP 23,614 18,519 (78.4%) 10,905 (58.9%) 7,511 (40.6%) 94 (0.5%) Table 2 Number of nodes, edges, trees, sites, and mutations in total for the combined tree sequence and separately for each chromo- some. Chromosome Nodes Edges T rees Sites Mutations TOTAL 95,612 693,524 31,925 61,260 1,031,672 1 11,195 73,825 4,429 8,245 119,709 2 9,584 67,879 3,543 6,549 106,885 3 8,568 55,794 3,181 6,126 91,245 4 9,647 81,314 3,195 6,198 116,737 5 6,159 45,476 1,952 3,815 69,731 6 8,228 58,815 2,719 5,442 82,933 7 6,776 49,100 2,204 4,239 63,656 8 7,572 58,206 2,428 4,577 79,862 9 6,620 46,421 1,971 3,977 80,205 10 6,513 45,489 1,930 3,718 77,165 11 8,671 70,554 2,631 4,988 86,870 12 6,079 40,651 1,742 3,386 56,674 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 20, 2025. ; https://doi.org/10.1101/2025.01.14.633033doi: bioRxiv preprint Rebollo et al. 2025 13 Individual Individual — 1.5 — 1.0 — 0.5 — 0.0 — -0.5 Figure 1 Genealogical nearest neighbors (GNN) of indica individuals in red and japonica individuals in blue. The GNN matrix was hierarchically clustered using the average method, shows the inherent population structure, and clearly separated the indica and japonica rice subspecies. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 20, 2025. ; https://doi.org/10.1101/2025.01.14.633033doi: bioRxiv preprint 14 Genomic Prediction with Ancestral Recombination Graphs in rice A B Subspecies japonica indica Figure 2 Differential local tree structure at genome positions: (A) the Drought and Salt Tolerance (DST ) locus associated with pan- icle length only in japonica; and (B) a locus associated with number of panicle secondary branches and number of spikelets per panicle secondary branch in both indica and japonica. The tree in (A) shows a very clear and deep separation between indica and japonica while the tree in (B) shows segregation in both indica and japonica. Table 3 Correlation between predicted mean genetic values from model 2 and between predicted specific genetic values from model 4 for different relationship matrices and cross-validation scenarios. Presented is the mean correlation across 1000 replicates of each scenario, with standard deviations in parentheses. Scenario Relationship matrix Model 2 Model 4 SRM BRM SRM BRM CVI J→I J PRM 0.78 (0.05) 0.82 (0.05) 0.74 (0.05) 0.76 (0.05) SRM 0.98 (0.01) 0.98 (0.01) PRM 0.70 (0.08) 0.72 (0.07) 0.70 (0.06) 0.73 (0.05) SRM 0.97 (0.01) 0.97 (0.01) CVI→J PRM 0.10 (0.15) 0.13 (0.15) 0.16 (0.12) 0.20 (0.12) SRM 0.88 (0.05) 0.84 (0.07) CVJ→I PRM -0.16 (0.19) -0.04 (0.19) -0.08 (0.15) 0.04 (0.13) SRM 0.74 (0.14) 0.80 (0.15) CVI→I PRM 0.79 (0.05) 0.81 (0.04) 0.77 (0.05) 0.82 (0.04) SRM 0.99 (0.00) 0.97 (0.01) CVJ→J PRM 0.67 (0.08) 0.68 (0.07) 0.62 (0.07) 0.63 (0.07) SRM 0.97 (0.01) 0.96 (0.01) Presented is the correlation ( r) between the predicted mean genetic values and between the predicted specific genetic values. Model 2 comprises a uniform covariance structure for the genetic values while model 4 comprises a compound symmetric covariance structure. IND - indica; JAP - japonica; PRM - pedigree-based relationship matrix; SRM - site-based relationship matrix; BRM - branch-based relationship matrix. CVI J→I J trained on and predicted both indica and japonica, CVI→J trained on indica and predicted japonica, CVJ→I trained on japonica and predicted indica, CVI→I trained on and predicted indica, and CVJ→J trained on and predicted japonica. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 20, 2025. ; https://doi.org/10.1101/2025.01.14.633033doi: bioRxiv preprint Rebollo et al. 2025 15 B A C Alternative allele frequency Alternative allele frequency Alternative allele frequency in indica from ALL Alternative allele frequency in japonica from ALL indicas from ALL japonicas from ALL ALL 0.00 0.25 0.50 0.75 1.000.00 0.25 0.50 0.75 1.00 0 5000 10000 15000 0 5000 10000 15000 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 Figure 3 Alternative (non-ancestral) allele frequency spectrum: (A) alternative allele frequency of the subset of indica and japonica from the genomic dataset with all individuals (ALL); (B) alternative allele frequency of the ALL genomic dataset; and (C) scatter plot between the alternative allele frequency of indica and japonica from the ALL genomic dataset. The histogram in (a) shows a U-shape distribution for each subspecies, as expected for neutral diversity , with an excess of the lower alternative allele frequency for japonica. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 20, 2025. ; https://doi.org/10.1101/2025.01.14.633033doi: bioRxiv preprint 16 Genomic Prediction with Ancestral Recombination Graphs in rice A B Age of nodes Age of first mutation at each site C Age of mutations 0 2500 5000 7500 10000 12500 0.0000 0.0005 0.0010 0.0015 0.0020Density Density 0.000 0.001 0.002 0.003 0.004 0.005 Density 0.0002 0.0000 0.0001 0.0003 0 2500 5000 7500 10000 12500 0 2500 5000 7500 10000 12500 Chromosome 1 2 3 4 5 7 6 8 9 10 11 12 Figure 4 Tree sequence inference for each chromosome: (A) estimated age of nodes; (B) estimated age of mutations; and (C) esti- mated age of first mutation at each site. The ages in (A) represent the age of ancestors in the trees while the ages in (C) represent the age at which sites became polymorphic, i.e., the age of the oldest ancestor above which we inferred a mutation. All ages were heavily right-skewed toward the present, with all chromosomes presenting a similar pattern. Table 4 Predictive ability of models 2 and 4 for different relationship matrices and cross-validation scenarios. Presented is the mean predictive ability across 1000 replicates of each scenario, with standard deviations in parentheses. The best-performing combina- tions are presented in boldface. Scenario Relationship matrix Number of individuals Model 2 Model 4 T raining set Prediction set Mean genetic values Mean genetic values Specific genetic values IND JAP IND JAP IND JAP IND JAP IND JAP CVI J→I J PRM 262 264 65 66 0.70 (0.03) 0.40 (0.02) 0.67 (0.03) 0.40 (0.02) 0.55 (0.02) 0.36 (0.01) SRM 0.68 (0.01) 0.47 (0.01) 0.63 (0.02) 0.47 (0.01) 0.58 (0.01) 0.40 (0.01) BRM 0.70 (0.01) 0.48 (0.01) 0.68 (0.01) 0.48 (0.01) 0.60 (0.01) 0.42 (0.01) CVI→J PRM 262 66 0.15 (0.03) 0.08 (0.03) 0.03 (0.01) SRM 0.11 (0.03) 0.10 (0.03) -0.10 (0.02) BRM 0.10 (0.02) 0.10 (0.03) -0.14 (0.01) CVJ→I PRM 264 65 -0.33 (0.03) -0.35 (0.03) -0.21 (0.02) SRM 0.43 (0.03) 0.46 (0.03) 0.28 (0.03) BRM 0.20 (0.04) 0.32 (0.04) 0.17 (0.02) CVI→I PRM 262 65 0.70 (0.01) 0.68 (0.01) 0.54 (0.01) SRM 0.69 (0.01) 0.66 (0.02) 0.60 (0.01) BRM 0.71 (0.01) 0.70 (0.01) 0.58 (0.01) CVJ→J PRM 264 66 0.42 (0.02) 0.42 (0.02) 0.28 (0.01) SRM 0.49 (0.02) 0.49 (0.02) 0.34 (0.01) BRM 0.49 (0.02) 0.49 (0.02) 0.35 (0.01) Presented is the correlation ( r) between the predicted mean genetic values and mean phenotypic values across environments and between the predicted specific genetic values and mean phenotypic values within environments. Model 2 comprises a uniform covariance structure for the genetic values while model 4 comprises a compound symmetric covariance structure. IND - indica; JAP - japonica; PRM - pedigree-based relationship matrix; SRM - site-based relationship matrix; BRM - branch-based relationship matrix. CVI J→I J trained on and predicted both indica and japonica, CVI→J trained on indica and predicted japonica, CVJ→I trained on japonica and predicted indica, CVI→I trained on and predicted indica, and CVJ→J trained on and predicted japonica. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 20, 2025. ; https://doi.org/10.1101/2025.01.14.633033doi: bioRxiv preprint Rebollo et al. 2025 17 A B C D BRMSRM BRM BRM Value -0.01 0.00 0.01 0.02 Value -1 0 1 2 Subspecies indica japonica BRM Subspecies indica japonica Subspecies indica japonicaindica vs japonica BRMBRM SRM SRM 2 3 1 0 -1 -0.01 0.00 0.01 0.020.0200.0160.012 1.5 2.0 2.5 3.0 Figure 5 Comparison of the site-based (SRM) and branch-based (BRM) relationship matrices: (A) heatmap of the SRM; (B) heatmap of the BRM; (C) comparison of the diagonal elements of each matrix; and (D) comparison of the off-diagonal elements. The figure shows that the SRM and BRM revealed a very similar population structure with elements of the matrices being highly correlated, although on a different scale. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 20, 2025. ; https://doi.org/10.1101/2025.01.14.633033doi: bioRxiv preprint 18 Genomic Prediction with Ancestral Recombination Graphs in rice 1 2 3 4 5 6 7 8 9 10 11 12 Model 2 GRM Model 2 BRM Model 4 GRM Model 4 BRM 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 Chromosome −log10(p−value) A B C D S S Figure 6 Log-transformed p-values for the predicted mean SNP site effects: (A) model 2 with the site-based relationship matrix (SRM); (B) model 2 with the branch-based relationship matrix (BRM); (C) model 4 with the SRM; and (D) Model 4 with the BRM. Model 2 comprises a uniform covariance structure for the genetic values while model 4 comprises a compound symmetric covari- ance structure. The number of significant SNP site effects was 13 for model 2 and 14 for model 4 with the SRM and 10 for model 2 and 5 for model 4 with the BRM. No SNP site effect was significant across all models and RMs. Table 5 Correlation between predicted SNP site effects from models 2 and 4 with the site-based (SRM) and branch-based (BRM) relationship matrices. Model Relationship Matrix Model 2 Model 4 SRM BRM SRM BRM 2 SRM 0.94 0.96 0.90 BRM 0.94 0.91 0.97 4 SRM 0.96 0.91 0.94 BRM 0.90 0.97 0.94 Model 2 comprises a uniform covariance structure for the genetic values while model 4 comprises a compound symmetric covariance structure. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 20, 2025. ; https://doi.org/10.1101/2025.01.14.633033doi: bioRxiv preprint Rebollo et al. 2025 19 ALL INDJAP 8,319 27,5529,849 168 15,540 14,8155,278 61,260 23,614 50,854 Supplementary Figure 1 Number of SNP sites unique and shared across the genomic datasets (ALL, both subspecies to- gether; IND, indica only; JAP japonica only). The number of SNP sites unique to each dataset is given outside the circles, while the number of shared SNP sites is given inside the cir- cles. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 20, 2025. ; https://doi.org/10.1101/2025.01.14.633033doi: bioRxiv preprint 20 Genomic Prediction with Ancestral Recombination Graphs in rice Supplementary Table 1 Summary of the fixed effects, random effects, and residual for models 1-7. All models included a fixed over- all mean and environment main effects as well as heterogeneous block and residual variances across environments. The models differed in their random genetic terms, as described below. Model Name Fixed effects Random effects Mean Env Genetic Non-genetic Residual Geno Env:Geno Env:Block Env:Block:Plot 1 Diagonal x x Diag Diag Diag 2 Uniform x x Unif Diag Diag 3 Compound symmetry x x Unif ID Diag Diag 4 Main effects + diagonal x x Unif Diag Diag Diag 5 Factor ana- lytic 1 x x FA 1 Diag Diag 6 Factor ana- lytic 2 x x FA 2 Diag Diag 7 Factor ana- lytic 3 x x FA 3 Diag Diag Env - environment; Geno - genotype; Unif - uniform; Diag - diagonal; ID - identity; FA k - factor analytic of order k. Supplementary Table 2 Residual log-likelihood (LogLik), Akaike Information Criterion (AIC), and percentage of variance explained (ve) for models 1-7 with different relationship matrices. Model LogLik AIC v e PRM SRM BRM PRM SRM BRM PRM SRM BRM 1 -171,666.2 -171,527.4 -171,487.0 343,572.5 343,294.8 343,214.0 2 -172,530.8 -172,408.4 -172,399.4 345,223.6 344,978.7 344,960.8 3 -171,477.2 -171,324.4 -171,288.9 343,118.4 342,812.7 342,741.9 43.7% 54.8% 47.6% 4 -171,380.8 -171,240.6 -171,206.3 343,003.6 342,723.2 342,654.6 35.5% 50.7% 42.7% 5 -171,334.7 -171,202.6 -171,169.5 342,989.4 342,725.3 342,659.0 49.6% 55.2% 46.3% 6 -171,245.4 -171,145.3 -171,113.8 342,888.8 342,688.6 342,625.5 67.5% 68.9% 64.5% 7 -171,204.9 -171,145.3 -171,069.6 342,883.8 342,707.9 342,613.2 81.8% 78.3% 80.2% PRM - pedigree-based relationship matrix; SRM - site-based relationship matrix; BRM - branch-based relationship matrix. Supplementary Table 3 Number of haplotypes in the drought and salt tolerance (DST ) gene and in the locus associated with the number of panicle secondary branches (NSB) and the number of secondary spikelets per secondary branch (NSSB) for indica (IND) and japonica (JAP). For both loci, the haplotype was defined as the combination of the alleles at the two SNP sites in the region. The allelic states of each SNP , i.e., ancestral (ANC) or alternative (ALT), are ordered by genomic position and separated by a hyphen. Percentages for each subspecies are presented within parentheses. Locus DST NSB & NSSB IND JAP IND JAP ANC-ANC 30 (4%) 1084 (98%) 208 (27%) 889 (80%) ANC-ALT 145 (19%) 3 (0%) 542 (71%) 8 (1%) ALT-ANC 585 (77%) 23 (2%) 8 (1%) 213 (19%) ALT-ALT 2 (0%) 0 (0%) 4 (1%) 0 (0%) .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 20, 2025. ; https://doi.org/10.1101/2025.01.14.633033doi: bioRxiv preprint Rebollo et al. 2025 21 Supplementary Table 4 Local trees and their SNP covering the zoom-in region, bordered by two SNP sites found significant in the genome-wide association study (chromosome 6, positions 29476724 to 29561694), and the number of genotypes with alternative (i.e., non-ancestral) state in indica (IND) and japonica (JAP) and corresponding percentage for each subspecies within parentheses, along with the standardised SNP site effect ( zi). Local tree SNP IND JAP z i Tree 1 S6_29476724 8 (1%) 366 (33%) +3.44 S6_29476748 175 (23%) 66 (6%) +0.35 S6_29476763 530 (70%) 399 (36%) +3.52 S6_29476787 531 (70%) 396 (36%) -3.55 Tree 2 S6_29480408 535 (70%) 400 (36%) -2.78 tree 3 S6_29480471 530 (70%) 397 (36%) +2.58 S6_29480530 529 (69%) 397 (36%) +2.53 Tree 4 S6_29516992 520 (68%) 27 (2%) -0.70 Tree 5 S6_29517004 5 (1%) 357 (32%) +2.17 Tree 6 S6_29531546 753 (99%) 457 (41%) +3.03 Tree 7 S6_29557666 751 (99%) 453 (41%) +3.11 S6_29557756 738 (97%) 99 (9%) +0.55 S6_29557803 13 (2%) 352 (32%) +2.28 S6_29561680 751 (99%) 456 (41%) +3.20 Tree 8 S6_29561694 751 (99%) 409 (37%) +3.28 Supplementary Figure 2 Differential local tree structure at the ending SNP covering the genome positions in (A) the Drought and Salt Tolerance (DST) gene, associated with panicle length only in japonica and (B) a locus associated with number of panicle sec- ondary branches (NSB) and number of secondary spikelets per secondary branch (NSSB) in both indica and japonica. The tree in (A) shows a very clear and deep separation between indica and japonica while the tree in (B) shows segregation in both indica and japonica. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 20, 2025. ; https://doi.org/10.1101/2025.01.14.633033doi: bioRxiv preprint 22 Genomic Prediction with Ancestral Recombination Graphs in rice Supplementary Table 5 Haplotypes and their number of genotypes in indica (IND) and japonica (JAP) and corresponding percentage for each subspecies within parentheses, along with the standardised haplotype effect ( hj) for each local tree covering the zoom-in region, bordered by two SNP sites found significant in the genome-wide association study (chromosome 6, positions 29,476,724 to 29,561,694). Haplotypes were defined as the combination of the allelic state (0 for ancestral and 1 for alternative) of all SNP sites found in each tree, ordered by genomic position. Local tree Haplotype IND JAP h j Tree 1 Chr6: 29,476,724 - 29,476,787 0000 55 (7%) 645 (58%) 0.00 0001 2 (0%) 0 (0%) -3.55 0010 3 (1%) 1 (0%) +3.52 0011 519 (68%) 32 (3%) -0.03 0100 173 (23%) 66 (6%) +0.35 0101 2 (0%) 0 (0%) -3.04 1010 0 (0%) 2 (0%) +3.83 1011 8 (1%) 364 (33%) +1.21 Tree 2 Chr6: 29,480,408 0 227 (30%) 710 (64%) 0.00 1 535 (70%) 400 (36%) -2.78 Tree 3 Chr6: 29,480,471 - 29,480,530 00 232 (30%) 712 (64%) 0.00 01 0 (0%) 1 (0%) +2.53 10 1 (0%) 1 (0%) +2.58 11 529 (70%) 396 (36%) +2.56 Tree 4 Chr6: 29,516,992 0 242 (32%) 1083 (98%) 0.00 1 520 (68%) 27 (2%) -0.7 Tree 5 Chr6: 29,517,004 0 757 (99%) 753 (68%) 0.00 1 5 (1%) 357 (32%) +2.17 Tree 6 Chr6: 29,531,546 0 9 (1%) 653 (59%) 0.00 1 753 (99%) 457 (41%) +3.03 Tree 7 Chr6: 29,557,666 - 29,561,680 0000 11 (1%) 652 (59%) 0.00 0001 0 (0%) 3 (0%) +3.20 0011 0 (0%) 1 (0%) +5.77 0100 0 (0%) 1 (0%) +0.55 1000 0 (0%) 1 (0%) +3.11 1001 0 (0%) 3 (0%) +3.27 1011 13 (2%) 351 (32%) +5.56 1101 738 (97%) 98 (9%) +3.68 Tree 8 Chr6: 29,561,694 0 11 (1%) 701 (63%) 0.00 1 751 (99%) 409 (37%) +3.28 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 20, 2025. ; https://doi.org/10.1101/2025.01.14.633033doi: bioRxiv preprint Rebollo et al. 2025 23 Supplementary Table 6 The 10 most frequent haplotypes, their number of genotypes in indica (IND) and japonica (JAP) and corre- sponding percentage for each subspecies within parentheses, and their standardised effect ( hj) for the zoom-in region, bordered by two SNP sites found significant in the genome-wide association study (chromosome 6, positions 29,476,724 to 29,561,694). Haplo- types were defined as the combination of the allelic state (0 for ancestral and 1 for alternative) of all SNP sites found in the zoom-in region, ordered by genomic position. Haplotype IND JAP h j 000000000000000 6 (0.8%) 622 (56.0%) 0.00 001111110111011 499 (65.5%) 23 (2.1%) +2.95 101111101110111 4 (0.5%) 335 (30.2%) +3.93 010000000111011 167 (21.9%) 9 (0.8%) +3.80 000000000111011 49 (6.4%) 2 (0.2%) +3.79 010000000111010 0 (0.0%) 46 (4.1%) +3.78 001111100111011 4 (0.5%) 2 (0.2%) +3.44 001111110110111 6 (0.8%) 0 (0.0%) +3.14 101111100110111 2 (0.3%) 4 (0.4 %) +3.87 010000000000000 0 (0.0%) 5 (0.5 %) +0.35 Other 51 25 (3.3%) 62 (5.6%) - Supplementary Table 7 Predictive ability of models 2 and 4 with the branch-based relationship matrix (BRM) built with three ef- fective population sizes (Ne) for two cross-validation scenarios in indica (IND) and japonica (JAP). Values presented are the mean predictive ability across 100 replicates of each scenario, with standard deviations provided in parentheses. Results of the best- performing Ne for each scenario and population are presented in boldface. Scenario Ne Model 2 Model 4 Mean genetic values Mean genetic values Specific genetic values IND JAP IND JAP IND JAP CVI J→I J 23 0.68 (0.07) 0.41 (0.09) 0.69 (0.07) 0.41 (0.09) 0.58 (0.05) 0.39 (0.05) 1,500 0.71 (0.06) 0.48 (0.09) 0.69 (0.06) 0.48 (0.09) 0.57 (0.05) 0.42 (0.05) 150,000 0.71 (0.06) 0.48 (0.09) 0.68 (0.06) 0.48 (0.09) 0.60 (0.04) 0.41 (0.05) CVI→J 23 0.03 (0.11) 0.03 (0.11) -0.08 (0.07) 1,500 0.10 (0.12) 0.10 (0.12) -0.14 (0.07) 150,000 0.11 (0.12) 0.11 (0.12) -0.14 (0.07) Presented is the correlation ( r) between the predicted mean genetic values and mean phenotypic values across environments and between the predicted specific genetic values and mean phenotypic values within environments. Model 2 comprises a uniform covariance structure for the genetic values while model 4 comprises a compound symmetric covariance structure. CVI J→I J trained on and predicted both indica and japonica; CVI→J trained on indica and predicted japonica. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 20, 2025. ; https://doi.org/10.1101/2025.01.14.633033doi: bioRxiv preprint 24 Genomic Prediction with Ancestral Recombination Graphs in rice B A C Minor allele frequency Minor allele frequency Minor allele frequency in indica from ALL Minor allele frequency in japonica from ALL indicas from ALL japonicas from ALL ALL 0.00 0.25 0.50 0.75 1.000.00 0.25 0.50 0.75 1.00 0 5000 10000 15000 0 5000 10000 15000 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 20000 20000 0.00 0.25 0.50 0.75 1.00 Supplementary Figure 3 Minor allele frequency spectrum. (A) Minor allele frequency of the subset of indica and japonica from the ALL genomic dataset that includes all individuals from both subspecies, (B) minor allele frequency of the ALL genomic dataset, and (C) comparison between the minor allele frequency between indica and japonica from the ALL genomic dataset. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 20, 2025. ; https://doi.org/10.1101/2025.01.14.633033doi: bioRxiv preprint Rebollo et al. 2025 25 Ne = 23 Ne = 23 Ne = 23 Ne = 1,500 Ne = 1,500 Ne = 1,500 Ne = 150,000 Ne = 150,000 Ne = 150,000 Age of nodes Age of nodes Age of mutationsAge of mutationsAge of mutations Age of sites Age of sites Age of sites 0.000 50 100 1500 0.005 0.010 0.015 density A B C density density density density density density density 0.000 0.050 0.075 0.100 0.025 50 100 1500 50 100 1500 0.00 0.02 0.04 0.0000 0.0010 0.0015 0.0020 0.0005 0.000 0.003 0.004 0.001 0.002 0.00050 0.00000 0.00025 0.00075 25000 5000 7500 10000 12500 25000 5000 7500 10000 12500 25000 5000 7500 10000 12500 0.0010 0.0000 0.0005 0.0015 density 0.0002 0.0000 0.0001 0.0003 0.00010 0.00000 0.00005 0.00015 Age of nodes 0.000080.000060.000040.000020.00000 0.000080.000060.000040.000020.00000 0.000080.000060.000040.000020.00000 D E F G H I 0.005 0.06 Supplementary Figure 4 Estimated age of nodes (A, B, and C), mutations (D, E, and F), and sites (G, H, and I) per chromosome with three effective population sizes (Ne). The age of nodes represents the age of ancestors in the trees. The age of sites represents the age at which sites became polymorphic, which equals the age of the oldest ancestor above which we inferred a mutation. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 20, 2025. ; https://doi.org/10.1101/2025.01.14.633033doi: bioRxiv preprint 26 Genomic Prediction with Ancestral Recombination Graphs in rice H A B C D E F G H I BRM Ne = 23 BRM Ne = 1,500 BRM Ne = 150,000 BRM Ne = 1,500BRM Ne = 23 BRM Ne = 150,000 SRM SRM SRM SRM SRM SRM BRM Ne = 1,500BRM Ne = 23 BRM Ne = 150,000 Subspecies indica japonicaindica vs japonica 2 3 1 0 -1 -0.01 0.00 0.01 0.02 -0.01 0.00 0.01 0.02 -0.01 0.00 0.01 0.02 0.00 0.01 0.02 0.030.00 0.01 0.02 0.030.00 0.01 0.02 0.03 2 3 1 0 -1 2 3 1 0 -1 2.5 2.0 1.5 2.5 2.0 1.5 2.5 2.0 1.5 Supplementary Figure 5 Comparison of the site-based (SRM) and branch-based (BRM) relationship matrices with three effective population sizes (Ne). Heatmap of the BRM (A, B, and C), comparison of the diagonal (D, E, and F), and off-diagonal (G, H, and I) elements of both matrices. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 20, 2025. ; https://doi.org/10.1101/2025.01.14.633033doi: bioRxiv preprint Rebollo et al. 2025 27 1 2 3 4 5 6 7 8 9 10 11 12 Model 2 GRM Model 2 BRM Model 4 GRM Model 4 BRM 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 Chromosome Absolute standardized marker effect A B C D S S Supplementary Figure 6 Standardised best linear unbiased predictions of the SNP site effects obtained with the site-based relation- ship matrix (SRM) from models 2 (A) and 4 (B), and the branch-based relationship matrix (BRM) from models 2 (C) and 4 (D). .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 20, 2025. ; https://doi.org/10.1101/2025.01.14.633033doi: bioRxiv preprint