# Source paper archive - Source title: HomoeoGWAS enables experiment-wide homoeolog-group interaction mapping in allopolyploid crops - Author(s): Cheng Si, Yujiang Wu, Zhu Sun, Cai Yao, Wei Li, Haiwang Zhang, Lihui Wang, Xuemei Sun, Qiwen Zhong, Shipeng Yang - Source landing URL: https://www.biorxiv.org/content/10.64898/2026.11.06.764877v1 - PDF URL: https://www.biorxiv.org/content/10.64898/2026.11.06.764877v1.full.pdf - PDF SHA-256: 4feb28fefb87b6052671271985154a0956e00ff3efa66669c72e58a0752ded55 - Converter: pymupdf4llm 0.0.27 - Limits: Automated PDF text extraction; reading order, equations, tables and formatting may be incomplete or incorrect. Images are omitted, and no OCR was performed. Markup and links are neutralized for safe downloading. Consult the original PDF for the authoritative version. --- [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) # 1 HomoeoGWAS enables experiment-wide homoeolog-group 2 interaction mapping in allopolyploid crops 3 Cheng Si, Yujiang Wu, Zhu Sun, Cai Yao, Wei Li, Haiwang Zhang, Lihui Wang, Xuemei Sun, 4 Qiwen Zhong* and Shipeng Yang* 5 Key Laboratory for Protection and Genetic Improvement of Qinghai-Tibet Plateau Germplasm 6 Resources (Co-construction by Ministry and Province), Academy of Agriculture and Forestry 7 Sciences, Qinghai University, Xining 810016, China 8 9 *Corresponding authors: Shipeng Yang (yangshipeng@qhu.edu.cn); Qiwen Zhong 10 (13997135755@163.com) 11 **Abstract** 12 Allopolyploid crops carry homoeologous gene copies from divergent subgenomes, an 13 evolutionarily defined set of candidate interactions that association studies do not test directly. 14 HomoeoGWAS uses these copies to define the unit of interaction testing, the homoeolog group, 15 and analyses all eligible groups as one family with experiment-wide error calibration, combining 16 three complementary genotype encodings per edge with a structure-aware parametric bootstrap. 17 Family-wise error was controlled in 80-group tests in wheat, cotton and rapeseed under several 18 null models and in full-family tests of all four crops under a structure-aligned global null (379– 19 17,404 groups). The omnibus outperformed a SNP-by-SNP scan when interaction spanned 20 several genotype dimensions. HomoeoGWAS detected homoeolog-group interactions for heading 21 date in wheat, flowering time in rapeseed and fibre length in cotton; five groups at four loci 22 retained support after regional conditioning. Rapeseed homoeolog coupling showed a weak but 23 reproducible association with homoeologous-exchange frequency. Genome-wide interaction 24 analyses took minutes. 1 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 25 **Introduction** 26 Several major crops, including bread wheat, cotton, oilseed rape and peanut [1–4], are 27 allopolyploids that combine two or three divergent subgenomes in one nucleus, so that many 28 genes are present as homoeologous copies [5–7] . These copies are not independent: their expression 29 is coordinated and frequently biased towards one subgenome [8–10], their dosage balance can shape 30 phenotype [11], and in some species, such as _Brassica napus_, recurrent homoeologous exchanges 31 reshuffle copies between subgenomes [3,12–15] . Allopolyploid genomes therefore carry a biologically 32 defined set of candidate genetic interactions among homoeologous copies. 33 Standard association workflows do not test this set directly. Genome-wide association 34 studies test one marker or one gene at a time, and the mixed models used in polyploid crops 35 either treat the genome as a single kinship [16–18] or model allele dosage within a marker [19] . 36 Interaction scans in their general form test all marker pairs [20,21], a search space that reached 1.96 × 37 10¹¹ pairs for the 627,255 SNPs mapped to the 17,404 rapeseed homoeolog groups studied here, 38 with a multiplicity and computational cost that make it impractical as a routine analysis [22,23] ; 39 recent marginal tests instead aggregate the interactions of each variant with the rest of the 40 genome [24,25] . Subgenome-partitioned variance components with additive-by-additive terms 41 improve genomic prediction in wheat [26,27], but a variance component estimates an aggregate 42 contribution; it does not test which homoeologous genes interact. What is missing is a calibrated, 43 genome-wide association test whose unit is the homoeolog group itself. 44 Such a test must meet three requirements. It must be sensitive to interactions of different 45 form, because homoeologs can act through dosage, through shared haplotype structure or through 46 more complex genotype combinations, and a test built for one form can be blind to another. It 47 must control the experiment-wide error rate across thousands to tens of thousands of groups in 2 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 48 panels with strong population structure, in which residual variance itself often follows that 49 structure. And it must be fast and reproducible enough to be applied routinely to whole genomes. 50 Here we present HomoeoGWAS, which uses the evolutionary relationship between 51 homoeologous copies to define the unit of interaction testing and turns the homoeolog groups of 52 an allopolyploid genome into one testable family with experiment-wide error calibration. Each 53 homoeolog group of two, three or four copies is tested through its pairwise edges with three 54 complementary genotype encodings combined into one omnibus P value, and error is controlled 55 across the whole family by a parametric bootstrap minimum-P procedure with a structure-aware 56 variance model. A companion workflow partitions additive variance by subgenome. We show 57 that the interaction test is calibrated under the tested null models on the real genomes of four 58 crops at family sizes from 80 to 17,404 groups, that it is powerful across interaction architectures 59 and complements SNP-by-SNP scans, and that it detects homoeolog-group interactions for 60 heading date in wheat, flowering time in rapeseed and fibre length in cotton. The same per-copy 61 features show that homoeolog coupling across _B. napus_ accessions has a weak but reproducible 62 association with population-level homoeologous-exchange frequency. Complete interaction 63 analyses of these panels take minutes, and the software is available with an interface for AI 64 assistants. 65 **Results** 66 **HomoeoGWAS turns homoeolog structure into a testable interaction family** 67 Allopolyploid genomes carry their own map of candidate interactions: the homoeologous 68 copies of each gene. HomoeoGWAS uses this map to organize interaction testing (Fig. 1). We 69 applied it to hexaploid bread wheat (827 Watkins accessions; 2,143 A–B–D homoeolog triads 70 passing a genotype-only design gate), tetraploid upland cotton (1,245 accessions; 3,757 A–D 3 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 71 pairs), _Brassica napus_ rapeseed (926 accessions; 17,404 A–C pairs) and tetraploid peanut (230 72 accessions; 379 A–B pairs) (Fig. 1a). The pairings were locally collinear for 95.0–99.9% of edges 73 in wheat, cotton and rapeseed, against at most 17% under partners permuted within the same 74 chromosome pair, and for all nine units detected below; the small peanut family was too sparse 75 for this check to discriminate (Supplementary Table 1). The interaction workflow tests predefined 76 homoeolog groups (Fig. 1c). A group of two, three or four copies contributes one, three or six 77 pairwise edges; each edge is tested with three encodings of the two gene blocks (minor-allele 78 burden product, first-principal-component product and a kernel-Hadamard term), and the 79 encoding and edge P values are combined by the Cauchy combination into one group P value 80 (group-omniB). All groups of a declared experiment form one family, whose family-wise error 81 rate (FWER) is controlled by a parametric-bootstrap minimum-P procedure (bootstrap-minP) 82 drawn from the fitted null model; in the primary null (smooth_pc4), sample-specific variance 83 weights, estimated from squared whitened residuals as a smooth function of the top four genomic 84 principal components, are applied to the fitted null covariance. Encoding and edge P values 85 localize the evidence within a group. A companion additive workflow builds one genomic 86 relationship matrix per subgenome, estimates their variance components jointly by multi-kernel 87 restricted maximum likelihood (REML) and reports the proportion of phenotypic variance 88 explained (PVE) by each subgenome with a leave-one-chromosome-out (LOCO) marker scan 89 (Fig. 1b). Every run binds its inputs by SHA-256 and stops rather than returning partial results 90 when a preflight check fails; the software runs from the command line or Python, with an 91 optional interface for AI-assisted analysis (Fig. 1d). 4 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 92 **Experiment-wide error control is calibrated on four crop genomes** 93 We first simulated null phenotypes on the real genotypes of the wheat, rapeseed and cotton 94 panels, each with a family of 80 homoeolog groups, and analysed them with the complete 95 pipeline, including per-phenotype REML refitting and bootstrap-minP (B = 2,000, α = 0.05; Fig. 96 2a; Supplementary Fig. 1; Supplementary Data 1, sheets S2-3 to S2-7). Under a kinship 97 preserving Gaussian null the empirical FWER was 3.9% in wheat (49 of 1,250; Wilson 95% 98 confidence interval 3.0–5.1%) and 4.6% in cotton (57 of 1,250; 3.5–5.9%) with the earlier 99 homoscedastic null, and 4.5% in rapeseed (56 of 1,250; 3.5–5.8%) with smooth_pc4. Nulls 100 carrying an additive or a PC1 main effect gave 1.4–4.4% (500 phenotypes each). 101 When residual variance followed population structure, the homoscedastic bootstrap null 102 over-rejected in wheat (FWER 0.157) and cotton (0.065), although not in rapeseed (0.017; 2,000 103 phenotypes per cell). We therefore adopted smooth_pc4, which applies structure-aware sample 104 weights to the fitted null covariance, as the primary null and validated it on fresh phenotypes 105 (Methods). Across five noise models (Gaussian, Student-t, 1% contamination at 6 s.d., and PC1 106 aligned and non-PC heteroscedasticity) and three panels, smooth_pc4 held the FWER between 107 0.021 and 0.052, with every Wilson upper bound at or below 0.068 (Fig. 2a; Supplementary Fig. 108 1b). In a separate threshold-transfer diagnostic with fixed whitening and no bootstrap, learned 109 thresholds also transferred to independent held-out phenotypes (3.8–6.2% across five null layers 110 and all statistics, with an omitted-kernel diagnostic reported separately; Supplementary Fig. 1a). 111 Finally, we tested calibration at the scale of the real families. In prespecified experiments, 112 null phenotypes combining a polygenic component (generator weight 0.30) with residual 113 variance aligned with the first two genotype principal components (weight 0.70) were reanalysed 114 over each complete family with the delivered pipeline (smooth_pc4, B = 2,000; 640 phenotypes 115 per panel; Supplementary Table 2). Any-group rejection was 1.4% for the 2,143 wheat triads 5 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 116 (one-sided 95% upper bound 2.4%), 3.1% for the 17,404 rapeseed pairs (4.5%), 1.1% for the 117 3,757 cotton pairs (2.0%) and 3.1% for the 379 peanut pairs (4.5%), with no failed analysis. 118 Families of 500 and 2,000 groups gave 1.2–5.8% (Supplementary Fig. 2). Error control was 119 therefore maintained from 80 to 17,404 groups under the tested global-null models. 120 **The omnibus improves detection of distributed interactions and complements SNP-by-SNP** 121 **scans** 122 To estimate the power a user would obtain, we planted one pairwise interaction per 123 phenotype in a unit drawn uniformly from the complete family, on the same background, and 124 reanalysed each phenotype with the delivered pipeline (Fig. 2b; Supplementary Data 1, sheet S2 125 2). At a nominal PVE of 10%, defined on the untransformed simulated phenotype, power was at 126 least 0.84 in 11 of the 16 panel-by-architecture combinations (six architectures in wheat, five in 127 each two-copy panel; 200 phenotypes each). The hardest architectures were an alternating-sign 128 signal that cancels an unsigned burden (wheat 0.63, cotton 0.62, rapeseed 0.05) and a single 129 causal SNP pair in wheat and rapeseed (0.58 and 0.25), and power was lowest in rapeseed, whose 130 family is by far the largest. 131 The three encodings were complementary. Each single encoding failed badly on some 132 architecture, such as the burden under alternating signs, whereas group-omniB fell more than 133 0.02 below a single encoding in only 32 of 288 cell-by-encoding comparisons in an 80-group 134 conditional grid (largest deficit 0.17; Fig. 2d; Supplementary Fig. 3; Supplementary Data 1, sheet 135 S2-9). The omnibus thus reduces the need to know the interaction architecture in advance. 136 Group-omniB was also more powerful than a SNP-by-SNP scan when interaction was 137 distributed over a group. We compared it with a regenie SNP × SNP product-term scan 138 aggregated to the group by minimum P, on the same 320 phenotypes per cell in 80-group families 6 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 139 (5% PVE, one causal group), with both methods calibrated to the same empirical error under 140 Gaussian nulls (held-out FWER 4.8–5.7%; Fig. 2c; Supplementary Fig. 4b; Supplementary Data 141 1, sheet S2-10). Group-omniB detected significantly more in 10 cells, mostly where the 142 interaction spanned several genotype dimensions or edges (wheat multidimensional kernel, 0.82 143 versus 0.22); four cells were ties, and the SNP × SNP scan was more powerful in the two cells 144 where one SNP pair carried the interaction (rapeseed, 0.59 versus 0.89). The two approaches are 145 therefore complementary, with the group test favoured when interaction is distributed over a 146 homoeolog group. 147 **Detections localize a homoeolog group and its region** 148 We next asked what a significant group identifies (Supplementary Figs. 3 and 4a; 149 Supplementary Data 1, sheets S2-8 and S2-11). In the whole-family power simulations, rejections 150 of non-causal groups fell mostly on the chromosomes carrying the causal unit’s copies, consistent 151 with linkage, and occurred on other chromosomes in 2.1% of replicates; in cells with power of at 152 least 0.80 the causal unit ranked first in 66–98% of replicates. Two controls examined specificity 153 in the 80-group conditional grid at 10% PVE. Strong additive effects concentrated in one group, 154 without interaction, raised the family-wise false-rejection rate to 8.0%, 9.0% and 12.8% in wheat, 155 cotton and rapeseed when phenotypes were scored with a fixed whitening matrix on the 156 untransformed scale and without the bootstrap. These rejections fell almost entirely on groups 157 other than the one carrying the additive effects (0 of 42 rejected groups in wheat, 0 of 46 in 158 cotton and 2 of 73 in rapeseed). With per-phenotype REML on the transformed phenotype, as in 159 the delivered pipeline, a rescoring of 200 stored phenotypes per panel gave 6.0% in wheat (95% 160 interval 3.5–10.2%) and 8.5% in cotton (5.4–13.2%), against 9.5% and 10.5% with fixed 161 whitening (Supplementary Fig. 4c). A strong additive effect concentrated in one homoeolog 7 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 162 group can therefore inflate the family-wise error above the nominal level, whereas genome-wide 163 additive or PC1 main effects (1.4–4.4%, end-to-end calibration) and planted additive variants of 164 up to 20% PVE outside the tested gene blocks (1.1–5.2%, whole families; leakage control below) 165 did not; the conditions are compared in Methods. An interaction between one copy of a group and 166 a gene outside it led to any-group rejection in 12.2–23.2% of replicates. Under the mispaired 167 control, 22% and 20% of rejected groups in wheat and cotton were the group containing the 168 reference copy and 16% and 18% the group of the external partner gene. A significant group 169 therefore reports interaction evidence for the group and its genomic neighbourhood, and 170 identifying the interacting partner requires follow-up of the group. 171 **Subgenome-stratified scans recover known loci and partition phenotypic variance** 172 The additive workflow recovered known loci in each panel (Fig. 3a; Supplementary Table 3). 173 For cotton fibre length (n = 1,245; genomic inflation factor 1.11), the scan recapitulated all four 174 loci previously reported for the same accessions at genome-wide (Bonferroni) significance. For 175 wheat heading date (days to ear emergence; n = 827; environment-adjusted line means with 176 records beyond 6 s.d. of the trial mean set to missing), 3 of the 29 loci reported for the Watkins 177 panel had a genome-wide significant SNP within ±5 Mb, against 0.18 expected (P = 0.002; 178 genomic inflation factor 1.39), and Vrn-B1 and Vrn-D1 reached genome-wide significance in the 179 2006 trial. For rapeseed flowering time (n = 926; genomic inflation factor 1.12), 3 of 7 known 180 flowering loci were recovered at P < 10⁻⁵, against 0.60 expected (P = 0.016); BnaFLC.A10 and 181 the BnaFLC.A03b/FRI region approached but did not reach Bonferroni significance (P = 9.7 × 182 10⁻⁸ and 6.1 × 10⁻⁸; threshold 4.95 × 10⁻⁸). In peanut (n = 230), hundred-seed weight had two and 183 seed length one SNP at Bonferroni significance, and AhARF2-2 was recovered at P = 1.8 × 10⁻⁶; 184 GEMMA detected at least as many significant SNPs on the same data. 8 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 185 Subgenome partitioning ranked the subgenomes by the variance they explain. Across 14 186 wheat traits with 200 bootstrap fits each (multistart-audited for four traits), NIR grain starch was 187 the only trait that passed the prespecified stability gate (Fig. 3b; Supplementary Table 4). The D 188 subgenome explained 0.465 of starch variance (95% bootstrap interval 0.274–0.604) and ranked 189 first in 98.0% of bootstrap fits; its share of the genetic variance (0.79) lay at the 96.8th percentile 190 of random chromosome groupings. Point estimates for the application traits are given in 191 Supplementary Fig. 5. 192 **Homoeolog-interaction scans in four allopolyploids** 193 We applied the interaction workflow to each application trait with one frozen family and one 194 bootstrap-minP correction per experiment (smooth_pc4, B = 2,000; Fig. 4f; Supplementary Table 195 5). For wheat the phenotype was the environment-adjusted heading date, inverse-normal 196 transformed within the pipeline. Four of 2,143 wheat triads (global P = 0.003, the bootstrap 197 minP-adjusted P value of the most significant group), two of 17,404 rapeseed pairs forming one 198 locus (0.0075) and three of 3,757 cotton pairs for fibre length (0.0245) were significant at 5% 199 FWER (units 1–9 in Fig. 4c–e). Cotton fibre elongation (global P = 0.152) and, in the small 200 peanut panel (230 accessions), hundred-seed weight (0.99) and seed length (0.63) yielded no 201 significant group. The fibre-elongation scan had a genomic inflation factor of 1.25 and is flagged 202 for calibration review. 203 The detected groups showed no strong nearby additive associations and were not unusually 204 aligned with the leading axes of population structure. The copy-resolved additive signal of each 205 unit is shown in Fig. 4d; no unit had an additive SNP with P < 10⁻⁵ within ±1 Mb, and their 206 alignment with population structure fell within the matched background (10th–66th percentiles; 207 Supplementary Table 6). In a prespecified control, additive variants planted outside every tested 9 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 208 gene block and correlated with the detected interaction features left the experiment-wide rejection 209 rate within the 7.5% tolerance in all nine panel-by-strength cells (1.1–5.2%), and no detected unit 210 was withdrawn (Supplementary Table 7). With the formal features held fixed, each unit kept 211 nominal support (P < 0.05) in at least 19 of 20 accession-deletion folds, and in wheat two triads 212 kept support in all four leave-one-environment reanalyses and two in three of four (Fig. 4e; 213 Supplementary Table 8). 214 Conditioning on regional genotype principal components, a post-selection robustness 215 analysis, separated two tiers of evidence. Five groups at four loci retained their evidence: the 216 wheat chromosome 5 and 6 triads (units 1 and 2), the two rapeseed A02 × C02 pairs at one locus 217 (units 5 and 6) and the cotton A12 × D12 pair (unit 7). The wheat chromosome 4 triad (unit 3) 218 was uncertain, and the wheat chromosome 1 triad (unit 4) and the cotton A05 × D05 and A08 × 219 D08 pairs (units 8 and 9) were not robust; we treat these as candidate regional signals. With the 220 homoscedastic null, six further wheat triads were rejected; they are listed as a sensitivity tier and 221 not claimed, because that null over-rejects under structure-aligned heteroscedasticity in this panel 222 (Supplementary Tables 9 and 10). 223 Among the robust signals, unit 1 is the strongest wheat detection (adjusted P = 0.003) and 224 carries orthology and domain support for a ZnF-UBP ubiquitin-specific protease [28] (rice OsUBP2), 225 making it the lead for follow-up. Fig. 4b illustrates the allelic pattern the test targets (Fig. 4a) 226 with unit 4, a regional candidate and the only detected triad with at least 26 accessions in every 227 genotype class: in raw phenotype means not adjusted for relatedness or population structure, 228 accessions in the higher PC1-defined genotype class at both the 1B and 1D copies (n = 26) 229 headed about 4 days earlier than expected from the two single-copy classes (difference of 230 differences −4.1 days; bootstrap 95% interval −8.6 to −0.7). Prespecified three-copy conditional 231 models of the same trait found no experiment-wide A × B × D effect (Supplementary Table 11). 10 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 232 **Homoeolog coupling is associated with homoeologous exchange in rapeseed** 233 The per-copy features used by the interaction test also reveal how homoeologs co-vary 234 across accessions, independently of any trait. For each edge we defined homoeolog coupling 𝑟 ! 235 as the absolute correlation across accessions between the two copies’ gene PC1 features after 236 adjustment for ten genotype principal components; because each accession contributes one 237 product term, coupling decomposes by accession, and its null distribution re-pairs each copy with 238 a non-homoeologous copy on another chromosome (Fig. 5a). Rapeseed coupling was highly 239 reproducible: its full-panel reliability lower bound of 0.704 exceeded the reliability threshold of 240 0.138 (median half-panel correlation 0.567; Fig. 5f). 241 Rapeseed coupling rose with the frequency of homoeologous exchange (HE) reported for 242 260 _B. napus_ accessions [29] . Across the 3,038 of 17,404 edges that map one-to-one to that study’s 243 anchors, the association was weak but reproducible (Spearman ρ = 0.136; within-chromosome 244 circular-shift permutation P ≤ 5 × 10⁻⁴), with mean 𝑟 ! rising from 0.146 for edges without carriers 245 to 0.241 for edges with more than 26 (Fig. 5c,d; Supplementary Fig. 6). It persisted after post hoc 246 adjustment for SNP count, minor-allele frequency and heterozygosity (partial ρ = 0.133) and was 247 positive in 9 of 10 A chromosomes, consistent with a contribution of recurrent homoeologous 248 exchange to homoeolog coupling across accessions. Filtering heterozygous calls removed most of 249 the genome-wide coupling excess over re-pairing, whereas the association with exchange 250 frequency was unchanged after adjustment for per-edge heterozygosity (ρ = 0.137). 251 Accessions carrying reported exchanges contributed disproportionately to the aggregate 252 coupling statistic at strongly coupled edges. In the 123 accessions shared with that study 253 (identical names and sequencing runs), the mean per-accession contribution across the mapped 254 accession–edge pairs was higher for carriers than for non-carriers (0.307 versus 0.162), exceeding 255 a null that reassigns whole HE profiles within ecotypes (P ≤ 0.001; Fig. 5b,e). In post hoc strata 11 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 256 the carrier excess was concentrated at strongly coupled edges (0.176 at 𝑟 ! ≥0.3, P = 0.0033, the 257 smallest value attainable with 300 permutations; 0.036 at 0.1–0.3, P = 0.07; −0.011 below 0.1, P 258 = 0.77). At strongly coupled edges (𝑟 ! ≥0.3), therefore, the accessions that carry an exchange 259 contribute disproportionately to the coupling; this is a population-level pattern rather than a 260 classifier of individual carriers (area under the curve 0.519 against 0.518 for permuted profiles; 261 Supplementary Fig. 7), and HE calls and genotypes share sequencing data. 262 Wheat and cotton showed weaker coupling with lower edge-level reliability (full-panel lower 263 bounds 0.052 and 0.050 against 0.094 and 0.111; a kinship-whitened cotton estimator passed 264 narrowly; Fig. 5f; Supplementary Table 12). In cotton, coupling exceeded the re-pairing null 265 (mean 0.039 versus 0.034; P = 0.001, also after kinship whitening) and passed the prespecified 266 within-cultivar-group test marginally (one-sided P = 0.042); in a post hoc diagnostic it 267 disappeared when principal components were re-estimated within groups, so it varies with finer 268 population substructure (Supplementary Fig. 8). In wheat, the genome-wide excess was small 269 (0.00105; one-sided 95% upper bound 0.0021), and associations with introgression, copy number 270 and expression co-variation were not separable from variant density or absent (Supplementary 271 Fig. 6). These panels thus bound the aggregate coupling signal at the precision their edge-level 272 measurements allow. 273 **Subgenome-aware and experiment-wide analyses run in minutes** 274 We measured wall time and peak memory on two nodes, each with two 64-core AMD EPYC 275 7742 processors, as medians of three repeats (four for mvMAPIT), with declared thread budgets 276 for every tool (Fig. 6; Supplementary Table 13). On wheat with 83.3 million markers and 827 277 accessions at 64 threads, the HomoeoGWAS scan with A, B and D kernels and LOCO took 17 278 min, REGENIE 20 min and GEMMA, with one pooled kinship, 10.5 h (Fig. 6a); wall time grew 12 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 279 with the number of accessions up to 5,000 simulated genotypes at 10 million markers (Fig. 6b). 280 On the 10-million-marker wheat scan, three subgenome kernels with LOCO cost about 18% more 281 than one kernel with LOCO (Fig. 6d), and wall time fell 5.5-fold from 1 to 256 threads with little 282 gain beyond 32 (Fig. 6e). At 5,000 simulated genotypes HomoeoGWAS used more memory than 283 GEMMA (52 versus 6.5 GB; Fig. 6f). 284 With bootstrap-minP (B = 2,000) and 128 workers, the complete interaction analysis took 7.3 285 min for wheat (2,143 triads), 6.7 min for rapeseed (17,404 pairs) and 3.4 min for cotton (3,757 286 pairs) (Fig. 6g). Experiment-wide calibration added about 11% over B = 19 (Fig. 6d), and the 287 wheat analysis took 6.7–7.7 min across 8–256 workers (Fig. 6e). At 17,404 rapeseed groups and 288 64 declared threads, tools that test other hypotheses took 3.0 h (mvMAPIT, marginal epistasis, on 289 the node used for single-locus scans), 4.3 h (PLINK 1.9, all SNP pairs) and 6.2 h (REGENIE, 290 within-group SNP pairs) (Fig. 6c); these figures compare cost under each tool’s configuration, 291 and the highlighted REGENIE single-locus and mvMAPIT measurements include repeats flagged 292 for interference from other jobs. 293 The software is available as a command-line tool, a Python library and an interface for AI 294 assistants, whose evaluation on 30 frozen analysis tasks is reported in Supplementary Fig. 9 and 295 Supplementary Table 14. 296 **Discussion** 297 HomoeoGWAS defines the unit of genetic-interaction testing by evolutionary homology: the 298 group of homoeologous gene copies that allopolyploidy created. Treating all such groups in a 299 genome as one family turns homoeolog interaction into a genome-wide association analysis with 300 experiment-wide error calibration. HomoeoGWAS adapts established P-value combination and 301 resampling procedures [30,31] to this biologically defined unit, which makes genome-wide 13 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 302 homoeolog interaction mapping feasible while accounting for multiplicity across the complete 303 family of eligible groups; each group is evaluated with encodings that capture dosage-like, 304 haplotype-like and multidimensional interaction. This reorganisation has three consequences 305 shown here. First, the family becomes small enough to calibrate by resampling from the fitted 306 null model: calibration held under Gaussian, heavy-tailed, contaminated and structure-aligned 307 heteroscedastic nulls in wheat, cotton and rapeseed, and under a structure-aligned global null over 308 the complete families of all four crops, up to the 17,404 A–C pairs of rapeseed. Second, the 309 omnibus is powerful without prior knowledge of the interaction architecture, because its 310 encodings are complementary: each single encoding failed on some architecture that the 311 combination detected. Third, the analysis is fast: complete interaction analyses of the benchmark 312 panels took 3.4–7.3 min with 2,000 bootstrap replicates and 128 workers, whereas tools testing 313 all SNP pairs or marginal epistasis took hours for the same groups. 314 The comparison with a SNP-by-SNP scan shows when a group test is the right tool. In 80 315 group families with thresholds calibrated empirically under Gaussian nulls, group-omniB 316 detected substantially more when simulated interaction was spread over several genotype 317 dimensions or edges of a homoeolog group, and the pairwise scan was more powerful when a 318 single SNP pair carried the signal. The two approaches are complementary: a group test 319 aggregates evidence that a pairwise scan dilutes across many tests, and a pairwise scan resolves a 320 concentrated signal that a group test averages. A natural workflow is therefore a group scan 321 across the genome, followed by SNP-level analysis inside significant groups. 322 The crop applications show what the test finds in practice. Nine units at eight loci were 323 significant at 5% experiment-wide error for heading date in wheat, flowering time in rapeseed 324 and fibre length in cotton. No additive SNP reached P < 10⁻⁵ within 1 Mb of the detected units, 325 none of the nine units met the prespecified withdrawal criterion in the additive-leakage 14 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 326 simulations, and each kept nominal support in at least 19 of 20 accession-deletion reanalyses. 327 Conditioning on regional genotype structure separated five robust groups at four loci — the 328 wheat chromosome 5 and 6 triads, the rapeseed A02 × C02 locus and the cotton A12 × D12 pair 329 — from an uncertain wheat chromosome 4 triad and three regional candidates (the wheat 330 chromosome 1 triad and the cotton A05 × D05 and A08 × D08 pairs). The wheat chromosome 5 331 triad, the strongest detection, carries a candidate ZnF-UBP ubiquitin-specific protease supported 332 by orthology and protein-domain evidence, providing a target for future homoeolog-specific 333 functional validation. 334 Homoeolog coupling adds a genome-level view that interaction scans alone do not give. In _B._ 335 _napus_, coupling between A and C copies was highly reproducible, rose with the population 336 frequency of homoeologous exchange and, at strongly coupled edges (𝑟 ! ≥0.3), was carried 337 disproportionately by the accessions that carry an exchange. This association supports a 338 population-level relationship between homoeolog coupling and reported exchange frequency, 339 rather than individual-level prediction of carriers, and the strongly coupled edges are candidates 340 for studying how exchange shapes trait variation [29,32] . In wheat and cotton, coupling was weaker 341 and its edge-level measurement less precise, so these panels bound the aggregate signal rather 342 than resolving it locus by locus. 343 These results define how HomoeoGWAS detections are read. Calibration was demonstrated 344 under the null generators tested, including heteroscedasticity in which sample-specific variance 345 weights vary smoothly with the leading genomic principal components; it does not extend to a 346 strong additive effect concentrated in one homoeolog group, and strong control of the family 347 wise error rate when some groups carry true interactions was not established. A significant group 348 locates interaction evidence in a homoeolog group and its genomic neighbourhood, and 349 identifying the interacting partner and the causal genes calls for follow-up such as allele-specific 15 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 350 expression, near-isogenic lines or homoeolog-specific mutants and editing [33,34] . Independent 351 population replication will be important for assessing the transferability of these signals. The 352 framework applies to allopolyploids with divergent subgenomes and established homoeolog 353 maps; autopolyploids, mixed-ploidy panels and direct higher-order interaction coefficients are 354 outside its scope. 355 Allopolyploid crops are being resequenced at scale [15,35–39], and their homoeolog maps 356 improve with each new assembly. HomoeoGWAS provides a calibrated, fast and reproducible 357 way to ask of these resources a question that single-marker association cannot: which homoeolog 358 groups act together to shape a trait. 359 **Methods** 360 **Panels, genotypes and homoeolog groups** 361 We analysed four allopolyploid diversity panels. Wheat ( _Triticum aestivum_, AABBDD): 827 362 phenotyped accessions of the Watkins landrace collection (ENA PRJEB71453) [36], with variants 363 on IWGSC RefSeq v1.0 coordinates and v1.1 gene models [1] . Cotton ( _Gossypium hirsutum_, 364 AADD): 1,245 accessions of CottonGVD [37] on the Gh_CRI v1 assembly. Rapeseed ( _Brassica_ 365 _napus_, AACC): 926 accessions of BnaGVD [38] with flowering-time records [40] on the Darmor-bzh 366 v4.1 assembly [3] . Peanut ( _Arachis hypogaea_, AABB): 230 accessions of PeanutPan [39] on the 367 NDH108 assembly. Biallelic SNPs were assigned to subgenomes by chromosome of origin 368 (wheat A/B/D chromosomes; cotton A01–A13/D01–D13; rapeseed A01–A10/C01–C09; peanut 369 Chr01–10 = A, Chr11–20 = B) and stored as one PLINK file set per subgenome. Variant filtering 370 was panel-specific: wheat single-locus scans used 83,296,047 SNPs with minor-allele frequency 371 (MAF) ≥ 0.01 and call rate ≥ 0.90; cotton inputs were prefiltered at MAF ≥ 0.10 and scanned at 372 MAF ≥ 0.05 (1,105,457 SNPs); rapeseed inputs were biallelic SNPs with MAF ≥ 0.05 and 16 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 373 missingness ≤ 0.20 (1,009,546 SNPs); peanut inputs were 587,385 SNPs with MAF ≥ 0.05. The 374 interaction features applied a further MAF gate of 0.01 within the analysed samples. 375 Homoeolog groups contain one gene per subgenome. Wheat triads were taken from the high 376 confidence 1:1:1 triads of Ramírez-González et al. [8] . Of the 17,036 high-confidence triads, 12,003 377 had at least three SNPs with minor allele frequency ≥ 0.01 in every copy in the Watkins panel; 378 11,920 of these had an estimable design, and a genotype-only design gate that used no phenotype 379 retained the 2,143 triads with an effective sample size of at least 20 for the conditional target, a 380 maximum single-accession information fraction of at most 0.20, a target residual ratio of at least 381 0.10 and at least five SNPs in every copy. For cotton (3,757 A–D pairs), rapeseed (17,404 A–C 382 pairs) and peanut (379 A–B pairs), homoeologs were identified as reciprocal best DIAMOND [41] 383 hits between subgenome proteomes (homoeogwas prep-homoeologs), restricted to homoeologous 384 chromosome numbers in cotton and peanut and searched across all A–C chromosome pairs in 385 rapeseed; pairs required at least three SNPs per copy. SNPs were mapped to genes within ±2 kb 386 of the gene or mRNA model (±1 kb for peanut; homoeogwas prep-snps), and each mapping file is 387 bound to its genotype file by SHA-256. As a descriptive, genotype-free check of the pairings, we 388 measured local collinearity: genes were ordered by start coordinate, and for each tested edge we 389 counted how many of the ten nearest tested edges along the first copy (five on each side) have 390 partners on the same chromosome within ten positions of its partner in the ordered list of tested 391 partners. Two nulls were used (200 permutations each, seed 20261104): partners permuted within 392 the same chromosome pair, which keeps every partner on its chromosome and the number of 393 pairs per chromosome pair, and partners reassigned freely within the subgenome pair. Per-panel 394 results are given in Supplementary Table 1. 17 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 395 **Phenotypes** 396 Wheat heading date was recorded as days to ear emergence in four field trials (2006, 2010, 397 2014, 2020) [36] . Line values were the least-squares means of the model value ~ environment + line, 398 which removes trial effects from lines observed in different trial subsets. For single-locus scans, 399 records more than 6 s.d. from their trial mean were set to missing before adjustment (one record); 400 the interaction analysis and the PVE atlas used the adjusted means without this rule, after inverse 401 normal transformation. Cotton fibre length and elongation were the best linear unbiased 402 predictors distributed with CottonGVD [37] . Rapeseed flowering time was taken from Wu et al. [40] . 403 Peanut hundred-seed weight and seed length were taken from the PeanutPan study. The other 13 404 wheat traits of the PVE atlas were single-environment values or least-squares means as above. 405 **Additive workflow** 406 For each subgenome _s_, a genomic relationship matrix 𝐾 " = 𝑍 " 𝑍 "# / ∑ $%&' ! 2 𝑝 "% (1 −𝑝 "% ) was 407 computed from the 𝑚 " variants of the subgenome, with 𝑝 "% the allele frequency of variant 𝑗 and 408 𝑍 " the mean-imputed dosages centred by 2𝑝 "% 42, using variants with minor-allele frequency 409 (MAF) ≥ 0.01, and scaled to trace _n_ . The model " 𝑦= 𝑋𝛽+ 7 𝑢 410 " + 𝑒, 𝑢 " ∼𝒩(0, 𝜎 "( 𝐾 " ), 𝑒∼𝒩(0, 𝜎 )( 𝐼) (1) ( 𝐾 " ), 𝑒∼𝒩(0, 𝜎 ) + 𝑒, 𝑢 " ∼𝒩(0, 𝜎 " " 411 was fitted by restricted maximum likelihood (REML) [43] with L-BFGS-B from multiple 412 starting points (one deterministic, the others log-uniform; seed 2026; ten by default, five for the ( 413 PVE atlas), keeping the fit with the highest likelihood. The PVE of subgenome 𝑠 is 𝜎 " /(∑ " " 414 𝜎 " [(] " + 𝜎 )( ); components with PVE below 10⁻³ are flagged as boundary estimates. Uncertainty was 415 quantified by a parametric bootstrap (200 replicates drawn from the fitted model and refitted), 416 giving percentile intervals and the bootstrap frequency with which each subgenome ranks first. 18 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 417 Marker association used the one-degree-of-freedom generalized-least-squares statistic 𝜒 [(] = 𝑈 [(] /𝐼 418 with variance components fixed at the genome-wide fit; in the leave-one-chromosome-out 419 (LOCO) scan [44,45] the kernel of the subgenome carrying the tested chromosome is rebuilt without 420 that chromosome’s markers, without refitting the variance components. The genomic inflation 421 factor is median(𝜒 [(] )/0.4549 [46] . 422 **Homoeolog-group interaction test** 423 Each group contributes all pairwise edges between its copies (1, 3 or 6 edges for 2, 3 or 4 424 copies); no higher-order coefficient is fitted. For each copy, SNPs with MAF ≥ 0.01 in the 425 analysed samples were mean-imputed and standardized; for genes with more than 150 such SNPs, 426 separate deterministic subsamples of 150 SNPs were drawn for the burden and for the principal 427 component features, with seeds derived by SHA-256 from the feature seed, subgenome, gene and 428 feature type. An edge was tested when both copies retained at least three SNPs and its interaction 429 columns were estimable given the covariates and main effects; non-estimable encodings were 430 omitted from the combination, and groups without an estimable edge were excluded from the 431 family. Three encodings of the interaction between copies X and Y were formed: (i) the product 432 of the minor-allele burdens (the mean of standardized SNPs oriented to the minor allele), (ii) the 433 product of the first principal-component scores of the two gene blocks and (iii) up to nine cross 434 products of up to three unit-variance principal-component scores per copy, tested jointly as fixed 435 effect interaction columns (a low-rank form of the Hadamard product of the two linear gene 436 kernels). 437 The phenotype was inverse-normal transformed, Φ [*'] {(𝑟−0.5)/𝑛} for rank 𝑟. The null 438 model 𝑦= 𝐶𝛽+ ∑𝑢 " " + 𝑒 (𝐶, intercept) was fitted by REML from three starting points with one 439 relationship matrix per subgenome (standardized SNPs with MAF ≥ 0.01, 𝑍 " 𝑍 "# /𝑚 scaled to 19 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 440 trace 𝑛, eigenvalues floored at 10⁻⁶), and the whitening matrix 𝑊= 𝑉 [T] [ *'/(] was formed. Each 441 encoding was then tested by a whitened generalized-least-squares nested F-test of its interaction 442 columns, conditional on the covariates and on the main effects of both copies (burden, PC1 or top 443 PCs, respectively), with numerator degrees of freedom equal to the rank added by the interaction 444 columns and denominator degrees of freedom _n_ minus the rank of the full design. The three 445 encoding P values of an edge were combined by the Cauchy combination with equal weights [30,47], 446 𝑇= 𝑚 [*'] ∑ $%&' tan {(0.5 −𝑝 % )𝜋} and 𝑃= 0.5 −arctan(𝑇)/𝜋 for 𝑚 P values, and the edge P 447 values of a group were combined in the same way to give the group P value (group-omniB). 448 Encoding and edge P values are reported to localize the evidence and are not part of the discovery 449 family. 450 **Experiment-wide error control** 451 All groups of a declared analysis form one family. Null phenotypes 452 𝑦,∗ = 𝐶𝛽[ + 𝑉T '/( 𝑧,, 𝑧, ∼𝒩(0, 𝐼), 𝑏= 1, …, 𝐵 (2) 453 were drawn from the fitted null model, already on the transformed scale, and analysed with 454 the same frozen whitening matrix. For each replicate the minimum P value over the family was 455 recorded; the adjusted P value of group 𝑔 is /01 = 1 + #{𝑏: min ∗ ≤𝑝 . } 2 [𝑝] [2,] 456 𝑝 . = 𝐵+ 1 (3) 𝑝 . 457 and groups whose P value falls below the ⌊𝛼(𝐵+ 1)⌋-th smallest null minimum are rejected 458 at family-wise level α (single-step minimum-P) [31] . Replicates with a non-finite statistic contribute 459 a minimum of zero, which is conservative. We used B = 2,000 and α = 0.05 throughout. 460 In the smooth_pc4 null, the bootstrap residuals are rescaled by a smooth variance function of 461 population structure. The squared whitened residuals of the null fit were regressed on an intercept 20 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 462 and the top four eigenvectors of the mean subgenome relationship matrix by a gamma 463 generalized linear model with log link (iteratively reweighted least squares, at most 50 iterations); 464 the fitted variances 𝑤, normalized to mean one, define 465 𝑦,∗ = 𝐶𝛽[ + 𝑤 '/( ∘i𝑉T '/( 𝑧, j, Cov(𝑦,∗ ) = 𝐷 3'/( 𝑉T 𝐷 3'/(, 𝐷 3 = diag(𝑤) (4) ,∗ = 𝐶𝛽[ + 𝑤 '/( ∘i𝑉T '/( 𝑧, j, Cov(𝑦,∗ 3'/( 𝑉T 𝐷 3 𝑦, ,∗ ) = 𝐷 3 466 The weights therefore rescale the whole fitted covariance, including its relationship 467 components, sample by sample; because they are estimated from whitened residuals, they are a 468 working approximation to this structure rather than a likelihood estimate, and their adequacy was 469 assessed by simulation under two forms of heteroscedasticity (below). The observed statistics are 470 unchanged. The homoscedastic null (𝑤= 1) is available as a sensitivity analysis. smooth_pc4 471 was selected after the unmodified bootstrap and two wild bootstraps (Rademacher and Mammen 472 weights) failed internally frozen criteria in a pilot, and was then validated on fresh phenotypes 473 with criteria frozen before analysis. 474 **Calibration and power experiments** 475 **Families and generators.** Calibration and conditional power used, per panel, an 80-group 476 family selected by stratifying on chromosome and marker-count quartile with a hash-based order 477 that used no phenotype; larger families of 500 and 2,000 groups nest the 80 groups. Simulated 478 signals were pairwise copy-by-copy products, standardized: minor-allele burden product, PC1 479 product, multidimensional kernel (sum of products of the first ≤ 3 PCs), single SNP pair, 480 alternating-sign SNP sums and, for wheat triads, burden products on two edges. For a target PVE 481 𝑞, the signal 𝑠 was standardized, the background was centred, projected orthogonally to 𝑠 and 482 standardized (𝑒 ⟂ ), and 𝑦= r𝑞 𝑠+ r1 −𝑞 𝑒 ⟂, so that PVE is the exact finite-sample share of the 483 untransformed phenotype variance. The additive-only control adds per-copy PC1 main effects 484 without interaction; the mispaired control is the burden of one copy times a gene of another group. 21 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 485 **End-to-end calibration.** In wheat and cotton, 1,250 Gaussian null phenotypes with the 486 panel’s kinship covariance and 500 each with an added additive or PC1 main effect were 487 analysed by the complete pipeline (per-phenotype REML, B = 2,000) with the homoscedastic 488 null; rapeseed was analysed in the same way with smooth_pc4. Failures counted as rejections, 489 and the prespecified criterion was a Wilson upper bound ≤ 0.075. 490 **Variance-model validation.** In wheat, cotton and rapeseed, the homoscedastic and 491 smooth_pc4 bootstrap arms were applied to the same fresh phenotypes, generated independently 492 of the development pilot, under Gaussian, Student-t (5 d.f., unit variance), 1% contamination at 493 ±6 s.d., PC1-aligned heteroscedastic (the whole kinship-preserving draw scaled by a variance 494 multiplier 1 + 3 minmax(PC1), matching the form assumed by smooth_pc4) and non-PC 495 heteroscedastic nulls (2,000 phenotypes per layer, 500 for the last). The frozen criteria were a 496 Wilson upper bound ≤ 0.075 for smooth_pc4 in each core layer and a paired power loss ≤ 0.05 in 497 12 positive cells (observed maximum 0.022). 498 **Threshold transfer.** For each panel, statistic and null layer — Gaussian and additive-only 499 (core), Student-t, PC1-heteroscedastic and contaminated (stress), and an omitted-kernel 500 diagnostic — the rejection threshold (100th smallest family minimum P among 2,000 calibration 501 phenotypes) was applied to 2,000 held-out phenotypes, using the conditional test with a fixed 502 whitening matrix and no bootstrap (B = 0). The same banks provided the 80-group power grid 503 (500 replicates per cell; 1 or 4 causal groups; PVE 2%, 5%, 10%) and the comparisons of group 504 omniB with each of the three single encodings in 96 positive cells (288 comparisons), each 505 statistic at its own threshold. 506 **Whole-family calibration and power.** Null phenotypes 𝑦= 𝑢+ 𝑒 with 𝑢∼𝒩(0,0.30𝐾) 507 and 𝑒∼𝒩(0,0.70𝐷), where 𝐾 is the mean subgenome relationship matrix (mean diagonal one) 22 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 508 and 𝐷= diag(𝑑 4 /𝑑 [‾] ) with 𝑑 4 = exp{0.4 𝑧(PC1 4 ) + 0.4 𝑧(PC2 4 )} for the first two eigenvectors of 509 𝐾 and 𝑑 [‾] the mean of 𝑑 4 (heteroscedasticity in the residual component only, a form different from 510 the one assumed by smooth_pc4), were analysed over each complete family with the delivered 511 pipeline (INT, per-phenotype REML, smooth_pc4, B = 2,000; 640 phenotypes per panel). The 512 prespecified criterion was a one-sided 95% Clopper–Pearson upper bound ≤ 0.075. For power, in 513 wheat, cotton and rapeseed (peanut entered the whole-family calibration but not the power 514 experiment), one interaction signal was added on this background to a unit drawn uniformly 515 among units whose copies each retained at least three SNPs (all units of each family): six 516 architectures in wheat and five in each two-copy panel, at PVE 2%, 5% and 10% (48 cells × 200 517 phenotypes = 9,600 analyses). The background was projected orthogonally to the signal so that 518 PVE is the exact share of the untransformed phenotype variance before inverse-normal 519 transformation. The signal used all gated SNPs of a gene, whereas the engine scores at most 150; 520 genes with more than 150 mapped SNPs were 4.5% of tested gene copies in wheat and at most 521 0.4% in the other panels. Non-causal rejections were classified, in order, as on a chromosome 522 carrying any copy of the causal unit, on another chromosome with the same number, or elsewhere. 523 Family scaling reanalysed the same 500 Gaussian and PC1-heteroscedastic phenotypes per panel 524 in wheat, cotton and rapeseed at 500 and 2,000 groups (smooth_pc4). 525 **SNP-by-SNP comparator.** In the 80-group families of wheat, cotton and rapeseed, regenie 526 4.1.2 [48] was run per phenotype with whole-genome ridge regression and leave-one-chromosome 527 out predictions (step 1) and a test of the product of each within-group cross-copy SNP pair with 528 both main effects (step 2, –interaction-snp –rare-mac 0 –no-robust); the group score was the 529 minimum pair P value. Group-omniB used the observed group P value (fixed whitening, B = 0). 530 Both methods were thresholded at the 100th smallest family minimum of 2,000 Gaussian 23 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 531 calibration phenotypes and evaluated on 2,000 held-out phenotypes, on additive-only and 532 mispaired controls and on the same 320 phenotypes per positive cell; paired power differences 533 have Newcombe 95% intervals. regenie excludes samples missing at a tested SNP, whereas 534 HomoeoGWAS mean-imputes. 535 **Additive-effect conditions.** The analyses with additive effects and no interaction differ in 536 the additive effect (a genome-wide genotype or PC1 main effect, per-copy PC1 effects of one 537 control group, or a planted variant outside the tested blocks), in the phenotype scale, in whether 538 whitening is fixed or refitted per phenotype, in the threshold source and in whether the bootstrap 539 is used; Supplementary Table 15 lists the conditions and rates side by side. 540 **Location of control rejections.** In an exploratory analysis of the stored conditional-grid 541 results, each rejected group under the additive-only and mispaired controls was classified as the 542 control group, the group of the mispaired partner gene, another group on the same homoeologous 543 chromosome number, or elsewhere (one control group, 10% PVE, group-omniB). 544 **Refit diagnostic.** In an exploratory analysis in wheat and cotton, stored phenotypes were 545 rescored with the fixed whitener on the untransformed or transformed response and with per 546 phenotype REML on the transformed response (200 phenotypes per control, 500 Gaussian; each 547 variant with its own 1,000-phenotype calibration). 548 All confirmatory simulation designs were prespecified: code and parameters were frozen 549 with hashes in version control before the corresponding results were generated, and amendments 550 affected execution only. These protocols were not publicly registered. Seeds were derived from 551 fixed roots by SHA-256 or numpy SeedSequence. 24 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 552 **Application analyses** 553 Each application trait was analysed with one family and one bootstrap-minP correction 554 (group-omniB; configured at the edge level in rapeseed and at the group level otherwise, which 555 coincide for two-copy groups; burden cap 150, minimum three SNPs per copy, MAF ≥ 0.01, 556 three PCs; smooth_pc4; B = 2,000; bootstrap and feature seed 2026). 557 **Known loci.** Recovery was assessed in fixed windows around loci reported for the same 558 panels or around cloned genes: cotton fibre-length loci of Peng et al. [37] (±250 kb), 29 heading-date 559 loci of Cheng et al. [36] (±5 Mb; expectation and P value from 1,000 joint per-chromosome circular 560 shifts), seven rapeseed flowering loci (±100 kb; Poisson-binomial expectation for random 561 windows) and the peanut seed genes AhARF2-2 and AhCKX6 (±100 kb). We report genome 562 wide Bonferroni significance (wheat 6.0 × 10⁻¹⁰, cotton 4.52 × 10⁻⁸, rapeseed 4.95 × 10⁻⁸, peanut 563 8.51 × 10⁻⁸) and, for rapeseed and peanut, P < 10⁻⁵. Wheat was also scanned for the 2006 trial 564 alone, and peanut counts were compared with GEMMA 0.98.5 [18] on the same data. 565 **PVE atlas.** Fourteen inverse-normal-transformed wheat traits were fitted with A, B and D 566 kernels built from LD-pruned SNPs (plink2 –indep-pairwise 1000kb 1 0.2; MAF 0.01; about 567 90,000–95,000 SNPs per subgenome), five starting values for the observed fit and 200 bootstrap 568 replicates with one start each; four candidate traits were re-audited with five starts per replicate. 569 The prespecified gate required convergence with at least 80% successful bootstrap refits, a 570 maximum PVE change ≤ 0.10 across MAF thresholds (0.01, 0.05) and transformations, the same 571 top subgenome in models with and without a homoeolog kernel (the elementwise product 𝐾 5 ∘ 572 𝐾 6 ∘𝐾 7 ), a top-rank frequency ≥ 0.80 in both models and bootstrap contrasts with the other 573 subgenomes excluding zero. The D share of starch genetic variance was compared with 1,000 574 random groupings of the 21 chromosomes into three sets of seven. 25 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 575 **Robustness.** Accession deletion used 20 deterministic folds and, for wheat, leave-one 576 environment reanalyses with the environment-adjusted mean refitted from the remaining trials; 577 the interaction features were held fixed and nominal P values are reported. In the additive-leakage 578 control, a variant with MAF ≥ 0.05 within ±5 Mb of a detected unit, outside every tested gene 579 block and maximally correlated with the unit’s burden and PC1 edge products, was added with 580 PVE 0, 10% or 20% to the whole-family null generator (640 phenotypes per cell; criterion as for 581 the whole-family null); a detected unit was to be withdrawn if rejected in at least 7 of 640 582 replicates. Regional conditioning reran each family with the top five genotype PCs per copy 583 within ±1 Mb (sensitivity: ten PCs within ±5 Mb) as covariates; a unit was retained when the 584 Clopper–Pearson upper bound of its adjusted P was below 0.05 in both windows and not robust 585 when the lower bound exceeded 0.05 in either. Three-copy conditional models of wheat heading 586 date (stable A × B × D after all lower-order terms on the line mean; environment × A × B × D 587 after all lower-order environment interactions on 3,099 line × environment observations with a 588 line-block Rademacher wild bootstrap; α = 0.025 each) were taken from their frozen analyses; the 589 second was re-executed from its frozen launcher. Candidate genes were annotated by orthology 590 and protein domains (InterPro [49] ), and expression was taken from the wheat developmental atlas [8] . 591 **Unit checks.** For each significant unit we recorded the strongest additive SNP P value within 592 each copy (gene ± 20 kb) and within ±1 Mb, and the absolute correlation between each copy’s 593 PC1 feature and the top five genotype principal components, expressed as a percentile among all 594 units of the family and among the 50 units nearest in gated SNP count and mean MAF 595 (Mahalanobis distance). 26 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 596 **Homoeolog coupling** 597 Per-copy PC1 features were replayed from the formal interaction runs (reproducing the 598 formal P values to 10⁻¹¹) and residualized on ten genotype principal components computed from 599 20,000 random gated SNPs. Homoeolog coupling 𝑟 ! of an edge 𝑒 is the absolute Pearson 600 correlation |𝑟 ) | of the two residualized features. With 𝑧 8 and 𝑧 9 the centred, unit-norm feature 601 vectors over all 𝑛 accessions of the panel, the per-accession contribution 𝑠(𝑒, 𝑖) = 602 sign(𝑟 ) ) 𝑛 𝑧 8,4 𝑧 9,4 averages to 𝑟 ! . In a kinship-whitened sensitivity analysis the features were 603 premultiplied by 𝑉 [*'/(] of a genome-wide relationship model before correlation [50] . The null re 604 paired each copy with a copy of another unit on a different chromosome number (1,000 re 605 pairings). Reliability was estimated from 20 splits of the panel into halves along k-means clusters 606 (k = 40) of the principal components, as the median covariate-conditioned Spearman correlation 607 between half-panel edge values; its lower bound (5th percentile of 1,000 5-Mb block bootstraps) 608 was projected to full-panel scale by the Spearman–Brown formula and compared with the 609 minimum reliability needed for 80% power to detect ρ = 0.13 in an external edge-level test. 610 Rapeseed homoeologous exchange (HE) segments for 260 accessions were taken from Wang 611 et al. [29] . A copies were linked to the study’s anchor genes through one-to-one BnPIR 612 orthogroups [51], and the HE frequency of an edge is the number of accessions with a segment 613 covering the anchor’s C-gene midpoint; exchange direction (A replacing C or C replacing A) was 614 taken from each segment’s annotation in that study. Associations were tested by Spearman 615 correlation against a within-chromosome circular-shift null: edges were ordered by A-copy 616 position, the HE frequencies were circularly shifted within each A chromosome (2,000 shifts), 617 and the one-sided P value is (1 + number of shifted ρ ≥ observed)/(2,000 + 1), so that 5 × 10⁻⁴ is 618 its floor; covariate-adjusted partial correlations included SNP count, MAF and heterozygosity. In 27 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 619 rapeseed, filtering heterozygous calls removes most of the genome-wide coupling excess over re 620 pairing, so the association with HE frequency was also tested after adjustment for per-edge 621 heterozygosity. For the 123 accessions shared with that study (identical names and sequencing 622 runs), the carrier statistic is the difference in mean contribution between carrier and non-carrier 623 accession–edge pairs, with a null that reassigns whole HE profiles among accessions of the same 624 ecotype (1,000 permutations). In a post hoc analysis, edges were stratified by coupling strength 625 (𝑟 ! < 0.1, 0.1 ≤𝑟 ! < 0.3 and 𝑟 ! ≥0.3), and the carrier difference was recomputed within each 626 stratum against 300 permutations of the same null, with P = (1 + number of permuted differences 627 ≥ observed)/301. Intervals for mean coupling by HE-frequency bin are percentile intervals of 628 1,000 bootstrap resamples of edges within each bin; intervals for the mean contribution of 629 carriers and non-carriers are percentile intervals of 1,000 resamples of whole accessions, each 630 resampled accession bringing all of its accession–edge pairs. External wheat and cotton tests used 631 introgression and copy-number frequencies [36,52], pan-genome gene classes [53], homoeolog 632 expression co-variation [54,55] and _G. barbadense_ introgression and presence–absence variation [56], 633 each against spatial circular-shift nulls. Cotton within-group tests used the 1,170 accessions 634 assigned to cultivar groups Cul1 ( _n_ = 449) and Cul2 ( _n_ = 721) by Xu et al. [57] : coupling was 635 recomputed within each group with ten group-specific principal components, its excess over 636 1,000 density-matched re-pairings (partner drawn from the same SNP-count quintile) was 637 combined across groups with inverse-variance weights, and the one-sided P value used the same 638 permutation index in both groups. Each primary test and statistic was frozen in version control 639 before it was computed; post hoc analyses are labelled as such. The reliability threshold was set 640 after the genome-wide coupling summaries were available, before any split-half statistic was 641 computed, and the reliability assessment is therefore exploratory. 28 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 642 **Computational benchmark** 643 Timing used two nodes, each with two 64-core AMD EPYC 7742 processors (1,007 GB and 644 503 GB memory; single-locus scans and mvMAPIT on the first, the other interaction analyses on 645 the second), one cell at a time, with /usr/bin/time -v and three repeats (median and range; the first 646 repeat with a cold cache; four repeats for mvMAPIT; 24-h cap). Repeats overlapping other heavy 647 jobs are flagged and reported. Single-locus scans of wheat heading date (827 accessions; 0.1, 1, 648 10 and 83.3 million nested random markers; 64 threads) compared HomoeoGWAS (A, B and D 649 kernels with LOCO; streaming in chunks of 200,000 markers), GEMMA 0.98.5 (one kinship; -gk 650 1, -lmm 1) [18], REGENIE 4.1.2 (step 1 block size 1,000) [48] and PLINK 2.0.0-a.6.9 (–glm, no 651 kinship) [58] . Panels of 2,500 and 5,000 genotypes were simulated for timing by copying 652 chromosome segments (mean 20 Mb) from real lines, with a phenotype of heritability 0.5. 653 Interaction analyses used 128 HomoeoGWAS workers with one native thread each; PLINK 654 1.9.0-b.8 –epistasis (all SNP pairs of the selected groups), REGENIE within-group SNP-pair tests 655 and mvMAPIT 2.0.4 [59,60] used 64 threads; the PLINK 1.9 time includes writing about 0.95 TB of 656 output to network storage. B = 19 is the smallest bootstrap size the engine accepts and was used 657 to estimate the cost of calibration. 658 **AI-assistant evaluation** 659 Thirty frozen tasks (10 routine analyses, 8 input traps, 7 invalid requests and 5 interpretation 660 questions; data cut from the real panels) were given to Claude Code CLI 2.1.296 with model 661 claude-sonnet-5-5 and Codex CLI 0.162.1 with model gpt-6-sol in isolated sandboxes with (A) 662 the HomoeoGWAS skill and tool server, (B) the command-line interface and documentation, or 663 (C) no tools (invalid-request and interpretation tasks only), three repeats each (432 runs). Routine 664 tasks were scored deterministically against expert reference outputs; other tasks combined 29 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 665 deterministic checks with an isolated claude-sonnet-5-5 judge, cross-checked on 311 runs by a 666 gpt-6-sol judge run through Codex CLI (agreement 0.913, κ = 0.818). Task definitions, prompts 667 and scoring rubrics are provided in the reproducibility archive. A before–after comparison of the 668 software without and with the added input guards used three input-trap tasks (36 runs per release) 669 and two valid control tasks (24 runs per release) across both agents and arms A and B (120 runs). 670 **Software and reproducibility** 671 HomoeoGWAS is written in Python and released under the MIT licence (GitHub, PyPI and 672 documentation site). A preflight refuses to run on missing genotype files, a SNP-to-gene map 673 whose genotype hash does not match, duplicated sample or gene identifiers, mismatched sample 674 order, fewer than ten overlapping samples or an undeclared marker encoding. Every formal run 675 records a pre-run manifest with the configuration hash, the git commit and a clean-tree check; 676 checkpoint manifests bind the family, phenotype, samples, relationship matrices, features and 677 null fit by SHA-256. Parallel interaction analyses use forked workers with one native linear 678 algebra thread each, and experiment-wide decisions were identical across worker counts in a 679 wheat check. The exact software version used for each analysis is recorded in the Zenodo archive. 680 Configurations, frozen analysis protocols, figure source data and per-replicate outputs will be 681 deposited in a versioned Zenodo archive. 682 **Data availability** 683 The genotype and phenotype data analysed in this study are publicly available. Genotypes 684 and heading-date records of the Watkins wheat collection are available from the European 685 Nucleotide Archive under accession PRJEB71453 [36] . Genotypes and fibre-quality phenotypes of 686 the cotton panel are available from CottonGVD (https://db.cngb.org/cottonGVD/) [37] . Rapeseed 687 genotypes are available from BnaGVD (https://rapeseed.biocloud.net/) [38], and flowering-time 30 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 688 records were taken from Wu et al. [40] . Peanut genotypes are available from PeanutPan 689 (https://cgm.sjtu.edu.cn/PeanutPan/), and seed-trait phenotypes were taken from Zhao et al. [39] . 690 Homoeologous-exchange segments of _B. napus_ were taken from Wang et al. [29] ; the additional 691 wheat and cotton datasets used for external comparisons are cited in Methods. The source data 692 underlying the figures, together with simulation configurations, analysis protocols and per 693 replicate outputs, will be deposited in a versioned Zenodo archive. 694 **Code availability** 695 HomoeoGWAS is open-source software under the MIT licence. The current release (v2.2.0; 696 Python ≥ 3.10) is available from GitHub (https://github.com/Shipeng-Yang/HomoeoGWAS, 697 release v2.2.0) and PyPI (https://pypi.org/project/homoeogwas/2.2.0/), with documentation at 698 https://shipeng-yang.github.io/HomoeoGWAS/. All releases are archived at Zenodo 699 (https://doi.org/10.5281/zenodo.20668367). Scripts that reproduce the simulations, application 700 analyses, AI-assistant evaluation and every figure will be made available at 701 https://github.com/Shipeng-Yang/HomoeoGWAS-reproducibility and archived in the same 702 Zenodo archive. 703 **References** 704 1. International Wheat Genome Sequencing Consortium (IWGSC). Shifting the limits in 705 wheat research and breeding using a fully annotated reference genome. _Science_ **361**, 706 eaar7191 (2018). 707 2. Hu, Y. et al. _Gossypium barbadense_ and _Gossypium hirsutum_ genomes provide insights 708 into the origin and evolution of allotetraploid cotton. _Nat. Genet._ **51**, 739–748 (2019). 709 3. Chalhoub, B. et al. Early allopolyploid evolution in the post-Neolithic _Brassica napus_ 710 oilseed genome. _Science_ **345**, 950–953 (2014). 31 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 711 4. Bertioli, D. J. et al. The genome sequence of segmental allotetraploid peanut _Arachis_ 712 _hypogaea_ . _Nat. Genet._ **51**, 877–884 (2019). 713 5. Van de Peer, Y., Mizrachi, E. & Marchal, K. The evolutionary significance of polyploidy. 714 _Nat. Rev. Genet._ **18**, 411–424 (2017). 715 6. Cheng, F. et al. Gene retention, fractionation and subgenome differences in polyploid 716 plants. _Nat. Plants_ **4**, 258–268 (2018). 717 7. Weiss-Schneeweiss, H., Schneeweiss, G. M. & Wickett, N. J. Polyploid genome 718 evolution across land plants: progress and perspectives. _Curr. Opin. Plant Biol._ **90**, 719 102857 (2026). 720 8. Ramírez-González, R. H. et al. The transcriptional landscape of polyploid wheat. _Science_ 721 **361**, eaar6089 (2018). 722 9. Glombik, M. et al. Rapid reprogramming and stabilization of homoeolog expression bias 723 in hexaploid wheat biparental populations. _Genome Biol._ **26**, 147 (2025). 724 10. Dorussen, D., Knight, E., Simmonds, J. & Borrill, P. Homoeolog expression in polyploid 725 wheat mutants shows limited transcriptional compensation. _New Phytol._ **249**, 16–23 726 (2026). 727 11. Birchler, J. A. & Veitia, R. A. Gene balance hypothesis: connecting issues of dosage 728 sensitivity across biological disciplines. _Proc. Natl Acad. Sci. USA_ **109**, 14746–14753 729 (2012). 730 12. Gaeta, R. T. & Pires, J. C. Homoeologous recombination in allopolyploids: the polyploid 731 ratchet. _New Phytol._ **186**, 18–28 (2010). 732 13. Hurgobin, B. et al. Homoeologous exchange is a major cause of gene presence/absence 733 variation in the amphidiploid _Brassica napus_ . _Plant Biotechnol. J._ **16**, 1265–1274 (2018). 32 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 734 14. Lu, K. et al. Whole-genome resequencing reveals _Brassica napus_ origin and genetic loci 735 involved in its improvement. _Nat. Commun._ **10**, 1154 (2019). 736 15. Song, J.-M. et al. Eight high-quality genomes reveal pan-genome architecture and ecotype 737 differentiation of _Brassica napus_ . _Nat. Plants_ **6**, 34–45 (2020). 738 16. Yu, J. et al. A unified mixed-model method for association mapping that accounts for 739 multiple levels of relatedness. _Nat. Genet._ **38**, 203–208 (2006). 740 17. Kang, H. M. et al. Variance component model to account for sample structure in genome 741 wide association studies. _Nat. Genet._ **42**, 348–354 (2010). 742 18. Zhou, X. & Stephens, M. Genome-wide efficient mixed-model analysis for association 743 studies. _Nat. Genet._ **44**, 821–824 (2012). 744 19. Rosyara, U. R., De Jong, W. S., Douches, D. S. & Endelman, J. B. Software for genome 745 wide association studies in autopolyploids and its application to potato. _Plant Genome_ **9**, 746 plantgenome2015.08.0073 (2016). 747 20. Cordell, H. J. Detecting gene–gene interactions that underlie human diseases. _Nat._ 748 _Rev. Genet._ **10**, 392–404 (2009). 749 21. Mackay, T. F. C. Epistasis and quantitative traits: using model organisms to study gene– 750 gene interactions. _Nat. Rev. Genet._ **15**, 22–33 (2014). 751 22. Purcell, S. et al. PLINK: a tool set for whole-genome association and population-based 752 linkage analyses. _Am. J. Hum. Genet._ **81**, 559–575 (2007). 753 23. Wei, W.-H., Hemani, G. & Haley, C. S. Detecting epistasis in human complex traits. _Nat._ 754 _Rev. Genet._ **15**, 722–733 (2014). 755 24. Fu, B. et al. A biobank-scale test of marginal epistasis reveals genome-wide signals of 756 polygenic interaction effects. _Nat. Genet._ **57**, 3175–3184 (2025). 33 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 757 25. Stamp, J., Smith, S. P., Weinreich, D. & Crawford, L. Sparse modeling of interactions 758 enables fast detection of genome-wide epistasis in biobank-scale studies. _Am. J. Hum._ 759 _Genet._ **112**, 2198–2212 (2025). 760 26. Santantonio, N., Jannink, J.-L. & Sorrells, M. Prediction of subgenome additive and 761 interaction effects in allohexaploid wheat. _G3_ **9**, 685–698 (2019). 762 27. Santantonio, N., Jannink, J.-L. & Sorrells, M. A low resolution epistasis mapping 763 approach to identify chromosome arm interactions in allohexaploid wheat. _G3_ **9**, 675–684 764 (2019). 765 28. Komander, D., Clague, M. J. & Urbé, S. Breaking the chains: structure and function of the 766 deubiquitinases. _Nat. Rev. Mol. Cell Biol._ **10**, 550–563 (2009). 767 29. Wang, T., van Dijk, A. D. J., Zhao, R., Bonnema, G. & Wang, X. Contribution of 768 homoeologous exchange to domestication of polyploid _Brassica_ . _Genome Biol._ **25**, 231 769 (2024). 770 30. Liu, Y. et al. ACAT: a fast and powerful _p_ value combination method for rare-variant 771 analysis in sequencing studies. _Am. J. Hum. Genet._ **104**, 410–421 (2019). 772 31. Westfall, P. H. & Young, S. S. _Resampling-Based Multiple Testing: Examples and_ 773 _Methods for p-Value Adjustment_ (Wiley, 1993). 774 32. Stein, A. et al. Mapping of homoeologous chromosome exchanges influencing 775 quantitative trait variation in _Brassica napus_ . _Plant Biotechnol. J._ **15**, 1478–1489 (2017). 776 33. Krasileva, K. V. et al. Uncovering hidden variation in polyploid wheat. _Proc. Natl Acad._ 777 _Sci. USA_ **114**, E913–E921 (2017). 778 34. Wang, Y. et al. Simultaneous editing of three homoeoalleles in hexaploid bread wheat 779 confers heritable resistance to powdery mildew. _Nat. Biotechnol._ **32**, 947–951 (2014). 34 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 780 35. Walkowiak, S. et al. Multiple wheat genomes reveal global variation in modern breeding. 781 _Nature_ **588**, 277–283 (2020). 782 36. Cheng, S. et al. Harnessing landrace diversity empowers wheat breeding. _Nature_ **632**, 783 823–831 (2024). 784 37. Peng, Z. et al. CottonGVD: a comprehensive genomic variation database for cultivated 785 cottons. _Front. Plant Sci._ **12**, 803736 (2021). 786 38. Yan, T., Yao, Y., Wu, D. & Jiang, L. BnaGVD: a genomic variation database of rapeseed 787 ( _Brassica napus_ ). _Plant Cell Physiol._ **62**, 378–383 (2021). 788 39. Zhao, K. et al. Pangenome analysis reveals structural variation associated with seed size 789 and weight traits in peanut. _Nat. Genet._ **57**, 1250–1261 (2025). 790 40. Wu, D. et al. Whole-genome resequencing of a worldwide collection of rapeseed 791 accessions reveals the genetic basis of ecotype divergence. _Mol. Plant_ **12**, 30–43 (2019). 792 41. Buchfink, B., Xie, C. & Huson, D. H. Fast and sensitive protein alignment using 793 DIAMOND. _Nat. Methods_ **12**, 59–60 (2015). 794 42. VanRaden, P. M. Efficient methods to compute genomic predictions. _J. Dairy Sci._ **91**, 795 4414–4423 (2008). 796 43. Yang, J., Lee, S. H., Goddard, M. E. & Visscher, P. M. GCTA: a tool for genome-wide 797 complex trait analysis. _Am. J. Hum. Genet._ **88**, 76–82 (2011). 798 44. Lippert, C. et al. FaST linear mixed models for genome-wide association studies. _Nat._ 799 _Methods_ **8**, 833–835 (2011). 800 45. Yang, J., Zaitlen, N. A., Goddard, M. E., Visscher, P. M. & Price, A. L. Advantages and 801 pitfalls in the application of mixed-model association methods. _Nat. Genet._ **46**, 100–106 802 (2014). 35 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 803 46. Devlin, B. & Roeder, K. Genomic control for association studies. _Biometrics_ **55**, 997– 804 1004 (1999). 805 47. Liu, Y. & Xie, J. Cauchy combination test: a powerful test with analytic _p_ -value 806 calculation under arbitrary dependency structures. _J. Am. Stat. Assoc._ **115**, 393–402 807 (2020). 808 48. Mbatchou, J. et al. Computationally efficient whole-genome regression for quantitative 809 and binary traits. _Nat. Genet._ **53**, 1097–1103 (2021). 810 49. Paysan-Lafosse, T. et al. InterPro in 2022. _Nucleic Acids Res._ **51**, D418–D427 (2023). 811 50. Mangin, B. et al. Novel measures of linkage disequilibrium that correct the bias due to 812 population structure and relatedness. _Heredity_ **108**, 285–291 (2012). 813 51. Song, J. et al. BnPIR: _Brassica napus_ pan-genome information resource for 1689 814 accessions. _Plant Biotechnol. J._ **19**, 412–414 (2021). 815 52. Heuberger, M. et al. Analysis of a global wheat panel reveals a highly diverse 816 introgression landscape and provides evidence for inter-homoeologue chromosomal 817 recombination. _Theor. Appl. Genet._ **137**, 236 (2024). 818 53. White, B. et al. De novo annotation reveals transcriptomic complexity across the 819 hexaploid wheat pan-genome. _Nat. Commun._ **16**, 8538 (2025). 820 54. Wang, X. et al. HebQTLs reveal intra-subgenome regulation inducing unbalanced 821 expression and function among bread wheat homoeologs. _Genome Biol._ **26**, 218 (2025). 822 55. Zhang, Z. et al. Population-scale gene expression analysis reveals the contribution of 823 expression diversity to the modern wheat improvement. _Nat. Commun._ **16**, 11133 (2025). 824 56. Yang, Z. et al. Graph pan-genome illuminates evolutionary trajectories and agronomic 825 trait architecture in allotetraploid cotton. _Nat. Genet._ **58**, 218–229 (2026). 36 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 826 57. Xu, Y. et al. Genomic reconstruction of upland cotton domestication uncovers staged 827 selection, gene flow, and flowering-time adaptation. _Proc. Natl Acad. Sci. USA_ **123**, 828 e2601246123 (2026). 829 58. Chang, C. C. et al. Second-generation PLINK: rising to the challenge of larger and richer 830 datasets. _GigaScience_ **4**, 7 (2015). 831 59. Crawford, L., Zeng, P., Mukherjee, S. & Zhou, X. Detecting epistasis with the marginal 832 epistasis test in genetic mapping studies of quantitative traits. _PLoS Genet._ **13**, e1006869 833 (2017). 834 60. Stamp, J., DenAdel, A., Weinreich, D. & Crawford, L. Leveraging the genetic correlation 835 between traits improves the detection of epistasis in genome-wide association studies. _G3_ 836 **13**, jkad118 (2023). 837 **Acknowledgements** 838 We thank Claude Code (Anthropic) for assistance with software development, code and 839 manuscript review, and language editing; the authors checked all outputs and take full 840 responsibility for the content. 841 **Author contributions** 842 S.Y. conceived the study, developed HomoeoGWAS, performed the analyses and wrote the 843 manuscript. Q.Z. co-supervised the study and revised the manuscript. C.S., Y.W., Z.S. and C.Y. 844 contributed to data collection and curation. W.L., H.Z., L.W. and X.S. contributed to the 845 interpretation of the results. All authors read and approved the final manuscript. 846 **Competing interests** 847 The authors declare no competing interests. 37 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 848 **Figure legends** 849 38 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 850 851 **Fig. 1 | HomoeoGWAS workflow.** Additive (blue track) and interaction (red track) analyses 852 share the inputs and one output. **a**, Panels (subgenome composition; accessions; tested 853 homoeolog groups) and a homoeolog group, one copy per subgenome. 𝐺 ", genotypes of 854 subgenome 𝑠; 𝑦, phenotype. **b**, Additive model: subgenome relationship matrices 𝐾 ", multi 855 kernel REML, subgenome proportions of phenotypic variance explained (PVE) and a leave-one 856 chromosome-out (LOCO) scan (illustrative values). **c**, Interaction test. Groups of two, three or 857 four copies (4x, 6x, 8x) give one, three or six edges; each edge is tested with the burden product 39 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 858 (𝑏 5 ⋅𝑏 7 ), the PC1 product (𝑧 5 ⋅𝑧 7 ) and a low-rank kernel-Hadamard term (𝐾 5 ∘𝐾 7 ). Encoding P 859 values are combined by the Cauchy statistic 𝑇 into edge P values (𝑝 56, 𝑝 57, 𝑝 67 ) and the group P 860 value 𝑃 . . Experiment-wide control: variance function (smooth_pc4) of the whitened residuals 𝑟 4 861 on genomic eigenvectors 𝑣 4;, null phenotypes 𝑦 '∗ … 𝑦 6∗ and the distribution of family minima '∗ … 𝑦 6∗ on genomic eigenvectors 𝑣 4;, null phenotypes 𝑦 '∗ 862 min 2 𝑝 2∗ with the family-wise level α. **d**, Hash-bound inputs, preflight checks, report and access 863 routes. 864 40 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 865 866 **Fig. 2 | HomoeoGWAS controls family-wise error and gains power across three crop** 867 **genomes.** **a**, Family-wise error rate (FWER; α = 0.05) in 80-group families. End-to-end: 868 complete pipeline under Gaussian (1,250 phenotypes), additive and PC1 main-effect nulls (500 869 each); open, homoscedastic null; filled, smooth_pc4. Variance model: five noise models; 870 magenta, smooth_pc4; grey, homoscedastic null. Threshold transfer: fixed whitening, no 871 bootstrap, 2,000 held-out phenotypes per layer. Bars, Wilson 95% intervals; dashed line, 0.05. **b**, 872 Power at interaction PVE 2%, 5% and 10% (light to dark), one planted interaction per phenotype 873 in the complete family; 200 phenotypes per point. **c**, Power difference between group-omniB and 41 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 874 a regenie SNP × SNP scan (320 phenotypes per cell; PVE 5%; matched empirical error); bars, 875 Newcombe 95% intervals; filled, interval excludes 0. **d**, Power gain of group-omniB over each 876 single encoding, mean of six conditional cells (PVE 2%, 5% and 10% × one or four causal 877 groups); points beyond ±0.35 drawn at the axis limit; magenta, omniB higher. 878 42 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 879 880 **Fig. 3 | HomoeoGWAS subgenome-stratified scans and subgenome partitioning of** 881 **phenotypic variance.** **a**, Report cards for wheat heading date (n = 827), rapeseed flowering time 882 (n = 926), cotton fibre length (n = 1,245) and peanut seed traits (n = 230). Manhattan plots, 883 LOCO scans coloured by subgenome; red dashed line, Bonferroni threshold; grey dotted line, P = 884 10⁻⁵; red circles, Bonferroni-significant SNPs (peanut); ticks below the axis, known loci (black, 43 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 885 recovered; grey, not recovered). Bars at right, PVE of each subgenome and the residual. **b**, 886 Subgenome PVE of 14 wheat traits: grey points, 200 bootstrap fits; coloured point and bar, 887 estimate and 95% interval; right, subgenome with the largest PVE and its bootstrap probability; 888 bold, trait passing the multistart stability gate. 889 44 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 890 891 **Fig. 4 | Homoeolog-interaction scans in four allopolyploid crops.** **a**, Illustrative non 892 additive patterns of two copies; dashed line, additive expectation; magenta, double combination. 893 **b**, Wheat triad 4 (1B × 1D): raw (unadjusted) mean heading date of PC1-defined genotype 894 classes relative to –/– (two-means classes of the gene PC1 score, sign fixed so that the largest 895 magnitude score is positive; –, lower; +, higher); error bars, 95% t interval of the class mean; n, 896 accessions per class; dashed line, additive expectation. **c**, Interaction scans of wheat (2,143 triads), 897 rapeseed (17,404 pairs) and cotton (3,757 pairs). Each column is a genome window whose 45 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 898 stacked squares give the −log10 P of its strongest group (one square per unit; windows with 899 −log10 P ≥ 1.5 or a significant group); magenta, groups significant at 5% FWER, copies joined 900 by a bar; red dashed line, experiment-wide threshold; 1–9, significant units. **d**, Strongest additive 901 −log10 P of each copy (coloured by subgenome) and interaction −log10 P of the group (diamond); 902 dashed line, interaction threshold. **e**, Nominal P of each unit with the formal features fixed in 20 903 accession-deletion folds (grey) and leave-one-environment reanalyses (triangles, wheat); 904 diamond, full data; dashed line, P = 0.05; k/N, reanalyses with P < 0.05 (N = 20; 24 in wheat). **f**, 905 Global bootstrap-minP P (tile; grey, not significant), significant over tested units and interaction 906 scan inflation factor λ (red, λ > 1.15) for each experiment. Units 1–9: wheat chromosome 5, 6, 4 907 and 1 triads; rapeseed A02–C02 pairs; cotton A12–D12, A05–D05 and A08–D08 pairs. 908 46 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 909 910 **Fig. 5 | Homoeolog coupling is associated with homoeologous-exchange frequency in B.** 911 **napus.** **a**, Schematic of a homoeologous exchange (HE) and of coupling 𝑟 !, the absolute 912 correlation across accessions of the structure-adjusted gene-PC1 features of the A and C copies, 913 which decomposes into per-accession contributions. **b**, Three HE hotspots: on A03, A09 and A02, 914 the 2-Mb A-copy window with the highest published HE frequency [29] among windows with ≥ 8 915 mapped edges, shown within 0–4 Mb. Rows, the 123 accessions shared with Wang et al., ordered 47 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 916 by ecotype and number of exchanges. Left, HE calls (blue, A replacing C; pink, C replacing A); 917 right, mean per-accession contribution in 0.5-Mb bins (white ≤ 0 to magenta ≥ 1); bottom, 918 carriers per edge and mean contribution of carriers (blue) and non-carriers (grey). **c**, Coupling in 919 1-Mb windows (magenta, upwards) with the 99th percentile of the re-paired null (dashed) and HE 920 frequency (grey, downwards; accessions of 260), A and C subgenomes. **d**, Edge coupling by HE 921 frequency (0, 1–5, 6–13, 14–26 and > 26 carriers; 3,038 edges): grey, distribution; magenta, 922 mean and 95% edge-bootstrap interval. **e**, Mean contribution of carriers (blue) and non-carriers 923 (grey) by coupling stratum (post hoc strata; 95% accession-bootstrap intervals). **f**, Split-half 924 reliability of edge coupling in 20 splits (grey; magenta, median; filled, 𝑟 ! ; open, kinship 925 whitened); dashed line, half-panel reliability needed for 80% power at ρ = 0.13. 926 48 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 927 928 **Fig. 6 | Computational cost of subgenome-aware scans and experiment-wide** 929 **homoeolog-interaction analysis.** Points, median of three repeats (four for mvMAPIT); bars, 930 range; open symbols, repeats overlapping another heavy job. **a**, Single-locus scan of wheat (827 931 accessions; 64 threads) against marker number: HomoeoGWAS (A, B and D kernels with 932 LOCO), GEMMA (one kinship), REGENIE and PLINK2 (no kinship; computational lower 933 bound). **b**, The same against the number of accessions at 10 M markers (2,500 and 5,000 934 simulated, right of the dotted line). **c**, Complete interaction analysis (B = 2,000; 128 workers) 935 against the number of B. napus groups, with tools that test other hypotheses: PLINK 1.9 (all SNP 936 pairs), REGENIE G×G (within-group SNP pairs) and mvMAPIT (marginal epistasis). **d**, Cost of 937 subgenome kernels and LOCO (top) and of the number of bootstrap replicates B (bottom). **e**, 938 Thread scaling of the single-locus scan and the interaction analysis. **f**, Peak memory in three 49 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 939 workloads. **g**, Complete interaction analysis per species (groups × accessions). Tool settings, 940 Supplementary Table 13. 941 50 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 942 **Supplementary Information** 943 **HomoeoGWAS enables experiment-wide homoeolog-group interaction mapping in** 944 **allopolyploid crops** 945 Contents: Supplementary Figs. 1–9; Supplementary Tables 1–15 (provided as 946 Supplementary_Tables_1-15.xlsx); Supplementary Data 1 (SupplementaryData1.xlsx); 947 Supplementary References. 948 **Supplementary Figures** 949 51 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 950 951 **Supplementary Fig. 1 | Calibration of the interaction test in wheat, rapeseed and cotton.** 952 **a**, Threshold transfer: held-out FWER of the conditional test (fixed whitening, B = 0) for group 953 omniB (filled circles) and the burden (open circles), PC1 (squares) and kernel (triangles) 954 encodings under six null layers; thresholds from 2,000 calibration phenotypes per layer, 2,000 955 held-out phenotypes per cell. **b**, Variance model: FWER of bootstrap-minP (B = 2,000) under 956 five noise models; open, homoscedastic null; filled, smooth_pc4. **c**, Q–Q plots of the global P 957 under Gaussian (black), additive (dark grey) and structure (light grey) nulls; shading, pointwise 958 95% interval; KS P, Kolmogorov–Smirnov test of the Gaussian-null global P against uniform. 959 Bars, Wilson 95% intervals; vertical line, 0.05; grey band, 0.035–0.065. 960 52 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 961 962 **Supplementary Fig. 2 | Family-size scaling under the global null.** **a**, Any-group rejection 963 at 500 and 2,000 groups (smooth_pc4, B = 2,000) under Gaussian (filled) and PC1 964 heteroscedastic (open) nulls, the same 500 phenotypes per layer at both sizes; bars, Wilson 95% 965 intervals; grey band, 0.035–0.065. **b**, Phenotypes rejected at 500 groups only, at both sizes or at 966 2,000 groups only. 967 53 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 968 969 **Supplementary Fig. 3 | 80-group conditional power grid.** Power (shading) of group 970 omniB (O, orange outline) and the burden (B), PC1 (P) and kernel (K) encodings for each 971 architecture, PVE and number of causal groups; 500 replicates per cell; each statistic at its own 972 empirical threshold. Right, fraction of group-omniB replicates with one causal group in which a 973 non-causal group was rejected. 974 54 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 975 976 **Supplementary Fig. 4 | Control architectures, comparator controls and refit diagnostic.** 977 **a**, Any-group rejection (shading) under the additive-only and mispaired-partner controls for 978 group-omniB (O) and the burden (B), PC1 (P) and kernel (K) encodings; 500 replicates per cell. 979 **b**, Held-out FWER (n = 2,000) and additive-only and mispaired controls (PVE 10%, n = 320) for 980 group-omniB (filled) and the SNP × SNP scan (open) at matched error. **c**, Group-omniB rejection 981 in wheat and cotton with a fixed whitening matrix on the untransformed (open circles) or 982 transformed phenotype (squares) and with per-phenotype REML on the transformed phenotype 983 (filled); held-out n = 500, controls n = 200. Bars, Wilson 95% intervals; grey band, 0.035–0.065. 984 55 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 985 986 **Supplementary Fig. 5 | Subgenome partition of the application traits.** REML PVE of 987 each subgenome kernel (colours) and the residual (grey) for the scans of Fig. 3a. 988 56 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 989 990 **Supplementary Fig. 6 | External tests of homoeolog coupling.** **a**, Association of edge 991 coupling 𝑟 ! with external measures [1–7] (Spearman or partial ρ; rank-biserial for the non-core gene 992 class). Filled, prespecified primary or replication test; open, post hoc density-adjusted reanalysis; 993 colour, species; grey bars, 95% range of the within-chromosome circular-shift null; right, 994 permutation P. **b**, Wheat: partial ρ with expression co-variation for 𝑟 ! (filled) and kinship 995 whitened coupling (open); grey bars, 95% null range. 996 57 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 997 998 **Supplementary Fig. 7 | Per-accession HE detection by coupling contribution.** ROC 999 curves of the per-accession contribution for HE carrier status in the 123 shared _B. napus_ 1000 accessions [1], per edge × accession (left) and per 2-Mb window × accession (right). Magenta, 1001 observed; grey line and band, mean and 95% range of 200 within-ecotype reassignments; dotted 1002 line, chance. 1003 58 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 1004 1005 **Supplementary Fig. 8 | Cotton coupling and population structure (post hoc).** **a**, Coupling 1006 excess over a density-matched re-pairing null in the whole panel and in cultivar groups Cul1 and 1007 Cul2 [8], with ten principal components from the whole panel (grey) or re-estimated within the 1008 group (magenta). **b**, Excess with 3, 10 and 20 within-group principal components. Bars, ±1.96 s.d. 1009 of 200 re-pairings. 1010 59 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 1011 1012 **Supplementary Fig. 9 | AI-assistant interface evaluation.** **a**, Task success rate of Claude 1013 Code and Codex CLI on 30 frozen tasks with the HomoeoGWAS skill and MCP tools (magenta), 1014 the command-line interface and documentation (dark grey) or no tools (light grey); 432 runs. **b**, 1015 Success rate by task category. **c**, Pass rate of input-trap (magenta, n = 36) and valid control tasks 1016 (grey, n = 24) without (v2.1.1) and with (v2.1.2) the added input guards; open circle, strict 1017 scoring. Bars in a and b, Wilson 95% intervals. 1018 **Supplementary Tables** 1019 Supplementary Tables 1–15 are provided as one workbook (Supplementary_Tables_1 1020 15.xlsx), one sheet per table. 1021 **Supplementary Table 1 | Local collinearity of the tested homoeolog families and of the** 1022 **nine detected units.** Support, number of the 10 nearest family edges along the first copy whose 1023 partners lie on the same chromosome within ±10 family positions of the edge’s partner. Nulls 1024 (200 permutations each, seed 20261104): free, partners reassigned among all edges of the same 1025 subgenome pair; constrained, partners permuted within the same chromosome pair. Bottom, the 1026 nine detected units. 1027 **Supplementary Table 2 | Whole-family global-null calibration of the four interaction** 1028 **families.** Null phenotypes y = u + e, u ~ N(0, 0.30K), e ~ N(0, 0.70D), D = diag(d_i / mean d), 60 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 1029 d_i = exp{0.4 z(PC1_i) + 0.4 z(PC2_i)}; delivered pipeline (INT, per-phenotype REML, 1030 smooth_pc4, B = 2,000, α = 0.05); prespecified criterion: one-sided 95% Clopper–Pearson upper 1031 bound ≤ 0.075. 1032 **Supplementary Table 3 | Recovery of known loci by the subgenome-stratified scans.** a, 1033 cotton fibre length; b, rapeseed flowering time; c, wheat heading date; d, peanut seed traits. 1034 Windows and expectations as in Methods. 1035 **Supplementary Table 4 | Subgenome PVE atlas of 14 wheat traits.** REML with A, B and 1036 D kernels; 200 parametric bootstrap fits per trait; screen, screening fits; multistart, audited fits; 1037 gate_passed, prespecified stability gate (Methods). 1038 **Supplementary Table 5 | Homoscedastic versus smooth_pc4 bootstrap null for every** 1039 **interaction analysis.** B = 2,000, seed 2026, α = 0.05, one bootstrap-minP family per analysis; 1040 smooth_pc4 is the primary null. Observed P values are identical under the two nulls; only the 1041 threshold changes. Cells: global P / threshold / rejected groups. 1042 **Supplementary Table 6 | Checks of the nine detected units against nearby additive** 1043 **signal and population structure, and candidate genes.** a, unit summary; b, percentiles among 1044 all units (pct_B) and among 50 matched units (pct_M); c, allele-frequency differentiation; d, class 1045 tests; e, candidate genes of the four wheat triads (smooth_pc4). 1046 **Supplementary Table 7 | Planted additive-leakage control and regional conditioning.** 1047 Planted additive variant at 0, 10 and 20% PVE (cells C0, C10, C20; 640 phenotypes each) and 1048 regional conditioning of the nine detected units (±1 Mb with five PCs per copy; ±5 Mb with ten 1049 PCs per copy). 1050 **Supplementary Table 8 | Accession- and environment-deletion reanalyses.** Nine detected 1051 units with the formal features fixed: 20 deterministic accession folds (all panels) and leave-one 1052 environment reanalyses (wheat); nominal P values. 61 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 1053 **Supplementary Table 9 | Wheat triads significant only under the homoscedastic null** 1054 **(sensitivity tier).** Triads significant under the homoscedastic null; the six marked “homoscedastic 1055 only” form the sensitivity tier. 1056 **Supplementary Table 10 | Wheat matched edge family and per-direction analysis.** 1057 Rejected edges of the wheat matched edge family (smooth_pc4) and analyses of each subgenome 1058 direction (AB, AD, BD) with the environment-adjusted phenotype. 1059 **Supplementary Table 11 | Wheat three-copy conditional models.** Prespecified three-copy 1060 conditional models of wheat heading date (days to ear emergence) over 2,143 triads; α = 0.025 1061 each. 1062 **Supplementary Table 12 | External tests of homoeolog coupling.** Associations of 1063 homoeolog coupling with external measures; prespecified primary or replication tests and post 1064 hoc analyses are labelled. Null mean and s.d. of the within-chromosome circular-shift null where 1065 available. 1066 **Supplementary Table 13 | Computational benchmark: timing, memory, thread budgets** 1067 **and tool versions.** Median and range of repeats; peak resident memory; declared threads; 1068 interference flags. Tools: HomoeoGWAS, GEMMA 0.98.5, REGENIE 4.1.2, PLINK 2.0.0-a.6.9, 1069 PLINK 1.9.0-b.8, mvMAPIT 2.0.4. 1070 **Supplementary Table 14 | AI-assistant evaluation tasks and outcomes.** Success by agent 1071 and arm (A, skill and tools; B, command line and documentation; C, no tools), by task category, 1072 and release-guard re-test without and with the added input guards. Cross-judge agreement 0.913 1073 (κ = 0.818). 1074 **Supplementary Table 15 | Analyses with additive effects and no interaction: conditions** 1075 **and family-wise rejection rates.** Conditions and family-wise rejection rates (any group rejected) 1076 of the analyses with additive effects and no interaction. 62 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.06.764877; this version posted November 9, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.06.764877) (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 [available under a CC-BY 4.0 International license.](http://creativecommons.org/licenses/by/4.0/) 1077 **Supplementary Data** 1078 **Supplementary Data 1 | Calibration, power and comparator tables.** Sheets S2-1 to S2-11; 1079 the README sheet lists the content of each sheet. 1080 **Supplementary References** 1081 1. Wang, T., van Dijk, A. D. J., Zhao, R., Bonnema, G. & Wang, X. Contribution of 1082 homoeologous exchange to domestication of polyploid _Brassica_ . _Genome Biol._ **25**, 231 1083 (2024). 1084 2. Heuberger, M. et al. Analysis of a global wheat panel reveals a highly diverse 1085 introgression landscape and provides evidence for inter-homoeologue chromosomal 1086 recombination. _Theor. Appl. Genet._ **137**, 236 (2024). 1087 3. Cheng, S. et al. Harnessing landrace diversity empowers wheat breeding. _Nature_ **632**, 1088 823–831 (2024). 1089 4. White, B. et al. De novo annotation reveals transcriptomic complexity across the 1090 hexaploid wheat pan-genome. _Nat. Commun._ **16**, 8538 (2025). 1091 5. Wang, X. et al. HebQTLs reveal intra-subgenome regulation inducing unbalanced 1092 expression and function among bread wheat homoeologs. _Genome Biol._ **26**, 218 (2025). 1093 6. Zhang, Z. et al. Population-scale gene expression analysis reveals the contribution of 1094 expression diversity to the modern wheat improvement. _Nat. Commun._ **16**, 11133 (2025). 1095 7. Yang, Z. et al. Graph pan-genome illuminates evolutionary trajectories and agronomic 1096 trait architecture in allotetraploid cotton. _Nat. Genet._ **58**, 218–229 (2026). 1097 8. Xu, Y. et al. Genomic reconstruction of upland cotton domestication uncovers staged 1098 selection, gene flow, and flowering-time adaptation. _Proc. Natl Acad. Sci. USA_ **123**, 1099 e2601246123 (2026). 63