Transcriptome-wide association analyses reveal the impact of regulatory variants on rice panicle architecture and causal gene regulatory networks

Rice is a global staple, and panicle architecture is a major determinant of grain yield. While several key genes (e.g., OSH1, LAX1, Gn1a, OsSPL14/IPA1, DEP1, FZP) have been cloned from mutants or bi-parental populations, their regulatory networks and favorable natural alleles remain incompletely understood. GWAS in rice has mapped panicle trait QTLs, but the long LD decay and complex genetic architecture hinder causal gene identification, with few significant loci despite high heritability. Regulatory variation, including cis- and trans-effects on gene expression, is a key driver of phenotypic diversity. Inspired by the human GTEx framework, the authors are building a Rice GTEx resource. This study focuses on young panicles (1–2 mm) from 275 varieties to: (1) identify genes whose expression associates with panicle traits; (2) examine selection on regulatory variants; (3) develop a cis- and trans-expression component-based association strategy to pinpoint causal genes at small-effect loci and infer regulatory networks; and (4) functionally validate key regulators affecting spikelets per panicle (SPP).

Prior research identified genes shaping rice panicle development and yield, including OSH1, LAX1, Gn1a, OsSPL14/IPA1, DEP1, and FZP. Some harbor cis-variants altering expression and panicle traits (e.g., IPA1, FZP). Multiple rice GWAS have reported panicle trait QTLs, but with few genome-wide significant loci and many sub-threshold peaks due to extensive LD and polygenicity. Human GTEx and plant eQTL studies demonstrate the value of population transcriptomes to link genetic variation to expression and traits. In maize and other crops, eQTLs and TWAS have aided trait dissection, though confounding via LD can yield false positives. There is limited understanding in crops of how domestication/breeding selected cis-regulatory variants and how such selection influenced agronomic traits. This work builds on these gaps by integrating TWAS, GWAS, eQTL, and a new cis/trans expression component framework in rice young panicles.

Study design and materials: A diversity panel of 529 representative rice varieties was used across four field experiments (different locations/years). Panicle phenotypes measured included spikelets per panicle (SPP), number of primary branches (NPB), and panicle length (PL) from five plants per variety. For transcriptomics, 300 varieties were selected for maximal genetic diversity; 275 yielded young panicle (1–2 mm) RNA-seq data. Zhenshan 97 was replicated across four plantings to assess batch effects. Phenotyping and integration: BLUP was used to integrate phenotypes across environments (year as random effect). Phenotypes were rank-normalized. Genotyping and GWAS: 10,567,425 variants (MAF>0.05) from RiceVarMap V2 were analyzed by EMMAX linear mixed models; kinship from all common variants (PLINK). The genome-wide significance threshold was 2.54×10⁻⁸ (≈1.97M independent tests) and suggestive threshold 5.07×10⁻⁷. SNP-heritability estimated via GCTA REML. The variance explained by loci at varying GWAS p-value thresholds (10⁻⁷ to 10⁻³) was quantified by LASSO (glmnet) with 10-fold CV R²; permutation tests assessed significance. RNA-seq and expression quantification: RNA from 1–2 mm young panicles was sequenced (>9 Gb PE150 per library). Reads were quantified with Salmon against RGAP gene models; miRNA primary transcripts were extended ±200 bp and quantified similarly. TPMs were log2(TPM+1) transformed; transcripts per gene were summed. PEER corrected the first three factors; expression per gene was rank-normalized. Expressed genes were those with expression >0.5 in >5% of varieties, yielding 30,869 genes (including miRNAs). TWAS: Linear mixed model (EMMAX) associated each gene’s expression with traits; FDR<0.05 defined significant associations. Gene set enrichment analyses (GSEApy) assessed enrichment for tissue/developmental expression programs, TFs, chromatin/epigenetic factors; TF family and motif enrichment via Plant Regulomics; GO enrichment performed. Linking GWAS loci to TWAS genes: For TWAS-significant genes, expression prediction from GWAS lead variants at varying p thresholds used LASSO; correlation (R²) between predicted expression and trait was computed. pQTL–eQTL hotspots were identified at relaxed thresholds (p<10⁻³ for both): a pQTL was a hotspot if it was an eQTL for many genes and those targets were enriched for >10 TWAS genes (BH-adjusted Fisher’s exact p<0.05). Selection on cis-regulatory variants: Derived alleles were defined using wild rice frequencies (RiceHap3). Lead cis-eQTLs (top variant per gene; p<10⁻⁵) within ±10 kb of TSS were analyzed for derived allele frequency (DAF) distributions compared to all variants. Effects of derived alleles on expression were stratified for positively/negatively associated genes (PAGs/NAGs) and by DAF bins; enrichment of GWAS signals at cis-eQTLs was assessed via Q–Q plots. Cis- and trans-expression component-based association study (cis-/trans-ECAS): For each gene, cis-expression component (cis-EC) was estimated via GCTA BLUP using variants within 100 kb to build the local GRM. Trans-expression component (trans-EC) was defined as residual expression after subtracting cis-EC. Traits were associated separately with cis-EC and trans-EC using LMMs; genes significant in both (cis-ECAS p<0.01; trans-ECAS FDR<0.01) were prioritized as putative causal. Similar ECAS was applied using expression of TWAS genes as molecular phenotypes (e-traits) to infer upstream regulators (cis-ECAS p<0.001; trans-ECAS FDR<0.001). Experimental validation: OsMADS17 function was tested via CRISPR/Cas9 knockout (ZH11 background). Panicle morphology and branching were quantified; statistics via two-tailed Student’s t-test. Dual-luciferase reporter assays tested OsMADS17’s transactivation on SDT (MIR156j) promoter fragments including a distal open chromatin region (dOCR, 4.2–4.4 kb upstream). Electrophoretic mobility shift assays (EMSA) evaluated binding. qRT-PCR quantified gene expression. Natural variation analysis identified OsMADS17 and SDT regulatory variants; haplotypes were related to expression and SPP. A breeding introgression introduced the favorable OsMADS17 haplotype (39-bp deletion) from IRAT109 into a temperate japonica line (1035); segregating progeny were phenotyped. Predictive modeling using favorable alleles: For cis-eQTLs of cis-/trans-ECAS genes, alleles up-regulating PAGs or down-regulating NAGs were labeled favorable. The number of favorable alleles per line was correlated with phenotypes in 254 varieties without transcriptomes; compared against random TWAS gene sets with significant cis-genetic variance.

- Panicle traits exhibited high SNP-heritability: SPP 0.852, NPB 0.820, PL 0.858. GWAS detected few genome-wide significant loci (at 2.54×10⁻⁸) but many sub-threshold peaks, indicating polygenic architecture with numerous small-effect loci. - TWAS in young panicles (275 varieties) identified 4,175, 5,844, and 6,839 genes significantly associated with SPP, NPB, and PL, respectively (FDR<0.05). Among top SPP-associated genes, multiple TFs (YABBY and MADS families) showed negative associations with SPP, consistent with roles in floral meristem determinacy. IPA1 (OsSPL14) expression positively correlated with SPP; miR156j (SDT) and miR529a were negatively associated with SPP and IPA1. - Enrichment analyses: SPP PAGs were enriched for genes preferentially expressed in young panicles and those down-regulated during panicle development; NAGs were enriched for genes up-regulated during development and for functions including flower development, meristem determinacy, transcriptional regulation, DNA/histone methylation, and miRNA biogenesis. TF family enrichments: PAGs (B3, SBP, NF-YA); NAGs (MADS, YABBY, TCP, C2H2). Promoter motif enrichments mirrored TF family results. - Expression of TWAS-significant genes is largely regulated by many small-effect loci: correlations between traits and LASSO-predicted expression (from GWAS pQTLs) increased markedly as p-value thresholds were relaxed (e.g., mean R² up to ~0.320 in the 275-panel). This supports a small-effect regulatory architecture. - At relaxed thresholds (p<10⁻³), numerous pQTL–eQTL hotspots were identified: 37 for SPP, 31 for NPB, and 176 for PL. Example SPP hotspot HSPP.05.1 was linked to 105 TWAS genes, including MIR529a, MOC1, and OsMADS1. - Selection on regulatory variants: Lead cis-eQTLs were enriched among high-DAF variants, indicating positive selection on cis-regulatory variants during domestication/breeding. For SPP, high-DAF derived alleles tended to up-regulate PAGs or down-regulate NAGs (i.e., favorable for SPP), while derived alleles with negative phenotypic effects had lower DAF. Correspondingly, cis-eQTLs with positive effects and high DAF showed stronger departures in GWAS Q–Q plots. - Cis-/trans-ECAS pinpointed putative causal genes despite weak single-variant GWAS signals. Identified 36 (SPP), 48 (NPB), and 99 (PL) cis-/trans-ECAS genes. Key examples: • SDT (MIR156j): both cis-EC and trans-EC negatively associated with SPP (p=4.39×10⁻¹⁰ and 1.06×10⁻¹⁰). Variants SDT-V1 (~5 kb upstream) and SDT-V2 (~24 kb upstream) jointly associated with SDT expression and SPP; haplotypes showed coordinated effects on expression and SPP. A ~2.7 kb intronic insertion in LD was noted. Despite modest GWAS p for SPP (3.53×10⁻⁵ for SDT-V1), integrating ECAS implicated SDT as causal. • OsMADS17 (AGL6): cis- and trans-EC negatively associated with SPP (p=1.61×10⁻³ and 5.46×10⁻¹⁰). Multiple cis-variants significantly affected expression, but local variants had weak direct SPP GWAS signals. - OsMADS17 regulates SDT: ECAS-based regulatory analysis predicted OsMADS17 as an upstream regulator of SDT. Transient dual-luciferase assays showed OsMADS17 transactivates SDT via a distal open chromatin region (dOCR, 4.2–4.4 kb upstream). Chromatin accessibility differences across haplotypes matched SDT expression. EMSA supported binding. Network analysis suggested additional negative regulatory links from OsMADS17 to RFL and OSH1. - Functional validation: CRISPR/Cas9 knockout of OsMADS17 increased SPP by 19.2% and number of secondary branches by 25.2%, with no significant change in NPB. - Natural variation and breeding utility: Two coding variants in OsMADS17 were not associated with SPP; 12 upstream cis-variants were, including a multi-allelic 11-bp/39-bp deletion (OsMADS17-V1) most strongly associated with expression (p=5.78×10⁻¹⁸). Deletions reduced OsMADS17 expression and increased SPP. Introgression of the 39-bp deletion haplotype from IRAT109 into a temperate japonica line increased NSB by 25.0% and SPP by 18.5% in segregating progeny. Favorable alleles were common in tropical japonica (OsMADS17 39-bp deletion freq 0.905) but rare in temperate japonica (0.087) and indica (0.028); SDT superior haplotype (G_G) similarly enriched in tropical japonica (0.687) and rare elsewhere. - Predictive power of cis-regulatory alleles: The count of favorable cis-eQTL alleles at cis-/trans-ECAS genes per variety correlated strongly with phenotypes in 254 lines lacking transcriptomes (e.g., SPP Pearson’s r=0.427, p=1.16×10⁻¹²), outperforming matched random TWAS gene sets. - Overall, cis-eQTLs of cis-/trans-ECAS genes showed stronger enrichment for GWAS signals than general TWAS genes, supporting their causal relevance.

The study addresses the challenge of dissecting polygenic panicle traits in rice by integrating GWAS, TWAS, eQTL, and a novel cis-/trans-expression component framework. Findings demonstrate that SPP and NPB are governed by numerous small-effect loci that largely act through gene regulatory variation. TWAS highlighted thousands of trait-associated genes, especially regulatory factors (TFs and epigenetic components) with coherent directional associations consistent with developmental biology. By decomposing expression into cis and trans components, the ECAS approach reduced LD-driven confounding and increased power to detect causal genes at small-effect loci where single-variant GWAS signals are weak. This uncovered 36 SPP causal candidates, including SDT (MIR156j) and OsMADS17, and enabled construction of causal regulatory links. Experimental assays confirmed that OsMADS17 up-regulates SDT via a distal regulatory element, forming a pathway whereby OsMADS17 increases SDT/miR156j, reducing IPA1 activity and lowering SPP. Selection analyses revealed that domestication and breeding favored cis-regulatory alleles that up-regulate positively associated genes or down-regulate negatively associated genes, effectively rewiring expression to increase grain number. The strong phenotype predictability from counts of favorable cis-eQTL alleles at ECAS genes underscores the practical breeding value of identified regulatory variants and networks.

This work provides a comprehensive regulatory view of rice panicle architecture by profiling young panicle transcriptomes, performing TWAS and GWAS, and introducing cis-/trans-ECAS to resolve causal genes and gene regulatory networks at small-effect loci. It identifies extensive regulatory variation shaping SPP/NPB/PL, reveals positive selection on beneficial cis-regulatory alleles, and validates a key regulatory axis (OsMADS17→SDT/miR156j→IPA1) affecting SPP. Superior cis-regulatory alleles at OsMADS17 and SDT are promising yet underused resources for yield improvement, as shown by CRISPR validation and introgression results. Future research should expand transcriptomic profiling across tissues and developmental stages, refine trans-EC estimation, functionally validate additional ECAS candidates, integrate chromatin and single-cell data to resolve direct targets, and apply the ECAS framework to other crops and complex traits.

- Association does not equal causation: TWAS and ECAS prioritize candidates, but many targets require further functional validation. - Trans-component estimation: Defining trans-EC as residual expression may include environmental/noise components and indirect effects. - Panel heterogeneity: Differences between the 529- and 275-variety panels suggest residual heterogeneity that can affect signal strength. - Tissue and stage specificity: Analyses focus on 1–2 mm young panicles; regulatory effects at other stages/tissues may be missed. - LD and mapping resolution: Long LD in rice can complicate fine-mapping of causal variants within hotspots. - Limited validation scope: Experimental validation centered on OsMADS17 and SDT; broader validation of network edges remains to be done.