Canalized gene regulatory networks stabilize floral polymorphism and enable modular transgressive expression

Canalization, Polymorphism, Stabilizing selection , Stellera chamaejasme , Transgressive , 20 Transcriptome 21 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint

22 Understanding how discrete phenotypic polymorphisms persist across populations despite ongoing 23 gene flow remains a central challenge in evolutionary biology and ecological genomics (Gray & 24 McKinnon, 2007; Kitano et al. , 2025) . In theory, recombination and migration should homogenize 25 phenotypic variation, yet many natural systems maintain distinct and stable morphs across 26 heterogeneous landscapes. Flower -color polymorphisms provide particularly powerful models for 27 addressing this ambiguity because they are visually discrete, ecologically consequential, and 28 repeatedly associated with pollinator behavior, abiotic stress, and reproductive isolation across diverse 29 plant lineages (Schemske & Bradshaw, 1999; Hopkins & Rausher, 2012; Wessinger et al. , 2014) . 30 Although phenotypic differentiation among floral morphs is well documented, the regulatory 31 mechanisms that allow alternative phenotypes to persist across spatially structured populations remain 32 poorly understood. 33 A key unresolved question is whether such polymorphisms are maintained primarily by selection 34 acting on individual loci or by stabilization of higher -order gene regulatory architectures. Eco -35 evolutionary theory predicts that when alternative phenotypes are adaptive, stabilizing or phenotype -36 associated selection may act on gene regulatory networks, producing canalized transcriptional states 37 with reduced variance around morph- specific optima (Waddington, 1942; Gibson & Wagner, 2000; 38 Masel & Siegal, 2009; Paaby & Rockman, 2014) . Under this framework, divergence among morphs 39 should be reflected not only in differences in mean expression, but also in reduced within-morph 40 variance and recurrent regulatory configurations across populations despite gene flow (Yeaman & 41 Whitlock, 2011; Tigano & Friesen, 2016; Ravinet et al. , 2017). Because complex traits are controlled 42 by coordinated interactions among many genes, such stability is expected to emerge at the level of 43 gene regulatory networks and co-expression modules rather than individual loci (Wagner, 2011; Boyle 44 et al., 2017; Payne & Wagner, 2019; Kadelka, 2026). 45 A major challenge, however, is distinguishing phenotype -associated regulatory structure from 46 patterns generated by neutral spatial processes. Geographic distance, demographic history, and 47 spatially autocorrelated environments can generate transcriptional divergence through drift or plasticity 48 alone, potentially leading to false inference of adaptive divergence if spatial structure is not explicitly 49 modeled (Meirmans, 2012; Wang & Bradburd, 2014; Rellstab et al. , 2015; Forester et al. , 2018; 50 Dauphin et al. , 2023) . Robust inference therefore requires analytical framework s that explicitly 51 partition transcriptomic variation into phenotype -associated (selection -consistent) and geography-52 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint associated (drift -consistent) components (Peres-Neto et al. , 2006; Capblancq et al. , 2020; Li et al. , 53 2025). If regulatory canalization underlies the persistence of floral morphs, transcriptomic similarity 54 should align more strongly with morph identity than with geographic proximity, producing a pattern 55 increasingly described as isolation by morph rather than isolation by distance . Perturbations of 56 stabilized regulatory systems provide a complementary approach for testing the limits of canalization 57 (Horta-Lacueva et al. , 2023; Runemark et al. , 2025) . When divergent regulatory backgrounds are 58 combined or reorganized, buffering may be relaxed, exposing epistatic interactions and generating 59 non-additive or transgressive expression (Dobzhansky, 1982; Rieseberg et al. , 1999; McManus et al. , 60 2010; Runemark et al. , 2025) . Importantly, such responses are often modular, with core regulatory 61 components remaining stable while specific subnetworks become destabilized (Gibson & Wagner, 62 2000; Wagner, 2008; Félix & Barkoulas, 2015; Kadelka & Murrugarra, 2024). 63 Floral pigmentation provides an ideal system for testing these ideas because it reflects integrated 64 metabolic and regulatory programs rather than the action of a single locus. Anthocyanin and carotenoid 65 biosynthesis pathways are tightly linked with primary metabolism, redox balance, stress responses, and 66 developmental regulation, making floral color a pleiotropic trait embedded within broader 67 physiological networks (Rausher, 2008; Albert et al., 2014; Smith, 2016) . Consequently, canalization 68 of pigmentation is expected to arise from stabilization of these broader regulatory systems, whereas 69 disruption of such networks may produce expression instability and phenotypes extending beyond the 70 parental range (Smith & Rausher, 2011; Wessinger & Rausher, 2015; Chen et al., 2025). 71 Here, we investigate these dynamics in Stellera chamaejasme L., a widespread alpine perennial of 72 the Qinghai–Tibetan Plateau and adjacent montane regions that exhibits four stable parental flower -73 color morphs of pure pink (PP), red–white (RW), yellow –white (YW), and pure yellow (YY) , 74 distributed across broad geographic and elevational gradients (Rana et al., 2024; Sun et al. , 2024). We 75 additionally identify a naturally occurring mosaic morph (RWh) in zones of sympatry between PP and 76 RW populations, providing a natural perturbation of regulatory architecture for testing whether gene 77 expression follows additive expectations or exhibits transgressive reorganization. Using a pan-78 transcriptomic framework integrating RNA -seq data across morphs, tissues, replicated populations, 79 and geographic regions, we test three predictions : (i) transcriptomic divergence is structured more 80 strongly by morph identity than by geography; (ii) parental morphs exhibit reduced within-population 81 variance consistent with transcriptional canalization; and (iii) the RWh morph reveals the limits of 82 regulatory stability through additive, dominant, or transgressive expression relative to parental morphs. 83 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint Together, this framework allows us to evaluate how canalized regulatory networks contribute to the 84 persistence of floral polymorphism and how its disruption generates transgressive regulatory states. 85

86 Study system for eco-evolutionary framework 87 Stellera chamaejasme is a widespread alpine perennial of the Qinghai –Tibetan Plateau and adjacent 88 montane regions. This species exhibits discrete flower -color polymorphism, including pure pink (PP), 89 red–white (RW), yellow –white (YW), and pure yellow (YY) (Rana et al., 2024). The persistence of 90 these morphs, despite potential gene flow, makes this system well suited for testing how stabilizing 91 selection structures regulatory architecture across populations. During field sampling, we identified a 92 naturally occurring red– white mosaic morph (RWh) in zones of sympatry between PP and RW 93 populations. These individuals frequently produced inflorescences bearing both RW -like and PP -like 94 floral elements, suggesting a probable hybrid origin from these two parental morphs. One hypothesis is 95 that RWh does not occupy an intermediate position in transcriptomic expression space but instead 96 exhibits extensive regulatory divergence, consistent with a strongly transgressive expression state. 97 Yellow-based morphs (YW and YY) were not observed in association with these mixed inflorescences 98 and were therefore treated as non -parental reference morphs in subsequent analyses. This system, 99 therefore, provides a natural eco -evolutionary framework to test whether floral morphs correspond to 100 stabilized transcriptional states, whether such states are canalized within populations, and how 101 hybridization and transgressive regulatory states may perturb stabilized regulatory architectures (Fig. 102 1). 103 Field sampling and phenotypic characterization 104 For each floral morph (PP, RW, YW, YY), three geographically distinct populations (at least 1 km far 105 apart) were sampled across Southwest China along elevational gradients (2,676–4,210 m a.s.l; Fig. 2A; 106 Supplementary Table S1). Within each population, flower tissue from three individuals was collected 107 to provide biological replication, with sampled plants spaced at least 100 m apart to minimize the 108 likelihood of sampling closely related or clonally derived individuals. The RWh morph was sampled 109 from two sympatric populations, reflecting its rarity. Floral tissue was collected at comparable 110 developmental stages to minimize ontogenetic effects on gene expression and was immediately flash -111 frozen in liquid nitrogen to ensure RNA integrity. Leaf, stem, and root tissues were additionally 112 collected to support transcriptome assembly. Geographic coordinates were recorded for all 113 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint populations, and voucher specimens were deposited at the Herbarium of the Kunming Institute of 114 Botany (KUN). 115 Floral spectral reflectance (300–700 nm) was measured for representative individuals using a fiber-116 optic spectrometer (USB2000+, Ocean Optics, Dunedin, FL, USA) equipped with a deuterium –117 halogen UV –VIS light source (DH -2000, Ocean Optics), with measurements calibrated against a 118 pressed barium sulfate (BaSO₄) pellet serving as a diffuse white reflectance standard. Spectral data 119 were processed in the R -package pavo2 v2.7.0 (Maia et al., 2019), providing quantitative phenotypic 120 validation of discrete floral morphs. 121 RNA sequencing, transcriptome assembly, and expression quantification 122 Total RNA was extracted and purified from the 1 mg of flower, leaf, stem, and root tissues using the 123 SpectrumTM Plant Total RNA Kit (STRN250, Sigma -Aldrich). RNA concentration was assessed using 124 the Qubit RNA Assay Kit with a Qubit 2.0 Fluorometer (Life Technologies, CA, USA) and purity was 125 assessed with a NanoDrop2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). 126 Strand-specific cDNA libraries were prepared using the NEBNext Ultra TM RNA/DNA Library Prep 127 Kit (New England Biolabs, USA) and sequenced on an Illumina HiSeq TM 2500 platform using paired-128 end (2 × 150 bp), generating raw reads of 19.58 to 33.66 million reads per library. In total, 160 129 libraries were generated (144 parental; 16 transgressive morphotype) (Table S2). Library preparation 130 and sequencing were performed by Novogene Bioinformatics Technology Co., Ltd. (Beijing, China). 131 A unified pan- transcriptome was constructed from all morphs and tissues. Raw paired- end reads 132 were quality filtered using Trimmomatic v0.39 (Bolger et al., 2014) to remove adapter sequences and 133 low-quality bases using sliding- window trimming (4 bp window, minimum Phred score of 20), with 134 reads shorter than 50 bp discarded. Following filtering, clean reads from flower, leaf, stem, and root 135 tissues were concatenated for each individual to generate a single paired -end dataset. These individual-136 level read sets were subsequently pooled across all samples for de novo assembly. The pan-137 transcriptome was assembled using Trinity v2.13.2 (Grabherr et al. , 2011; Haas et al. , 2013) with 138 strand-specific parameters. In addition to the global assembly, tissue -specific (flower, leaf, stem, and 139 root) and morph- specific (PP, RW, RWh, YW, and YY) assemblies were generated by pooling reads 140 within tissue or morph categories to evaluate assembly consistency and transcript representation across 141 biological groups. Redundant transcripts were removed using CD -HIT-EST v4.8.1 (Fu et al. , 2012) 142 using a 95% identity threshold to generate a non -redundant reference assembly. Assembly contiguity 143 was evaluated using Trinity’s expression- weighted ExN50 statistics, and completeness was assessed 144 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint with BUSCO v5.4.7 (Simão et al. , 2015) in transcriptome mode against the embryophyta_odb10 145 lineage dataset, reporting the proportions of complete (single -copy/duplicated), fragmented, and 146 missing conserved plant orthologs. 147 Coding regions were predicted from the non- redundant transcriptome using TransDecoder 148 v5.7.1(Haas et al. , 2013). Functional annotation was performed by similarity searches against the 149 UniProt/SwissProt protein database using BLASTX (Camacho et al., 2009) with an E-value threshold 150 of 1 × 10⁻⁵. Predicted protein sequences were further annotated using eggNOG-mapper v2 151 (Cantalapiedra et al., 2021) to assign Gene Ontology (GO) terms, KEGG pathways, and orthologous 152 functional categories based on precomputed orthologous groups. 153 A Salmon index was constructed from the non- redundant pan-transcriptome assembly. Transcript 154 abundance was estimated using Salmon v1.8.0 (Patro et al., 2017) in quasi-mapping (mapping-based) 155 mode, with reads mapped directly to the transcriptome without prior alignment. Transcript -level 156 abundance estimates were summarized as transcripts per million (TPM) and raw estimated counts. 157 Gene-level counts were generated using Trinity gene –transcript mappings. For differential expression 158 analyses, raw gene-level counts were normalized using the median -of-ratios method in the R -package 159 DESeq2 v1.42 (Love et al. , 2014) . Differential expression was assessed using generalized linear 160 models, with multiple -testing correction applied using the Benjamini –Hochberg false discovery rate 161 procedure. For co- expression and eco- evolutionary analyses, expression data were normalized using 162 trimmed mean of M -values (TMM) normalization in the R -package edgeR v3.44 (Robinson et al. , 163 2010) and log2-transformed. 164 All expression processing steps were performed using reproducible workflows implemented in R 165 v4.3, with key R -packages including DESeq2 v1.42, edgeR v3.44, limma v3.58, and tidyverse v2.0. 166 Normalization choices were tailored to the analytical objective of inferring gene -level inference for 167 differential expression and module -level inference for regulatory and eco -evolutionary analyses, 168 ensuring methodological consistency while maximizing biological interpretability. 169 Early test of transgressiveness in floral expression and sequence space 170 Field observations showed that the RWh exhibits a mosaic floral phenotype combining traits of the 171 sympatric morphs PP and RW, with several trait values extending beyond the parental range. Based on 172 this pattern, RWh was presumed to represent a transgressive hybrid, which we evaluated using 173 integrated orthologous phylogenomic and expression- based analyses. To assess the evolutionary 174 placement of the RWh morphotype relative to other morphs, a morph- level phylogenetic analysis was 175 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint conducted based on genome -wide orthologs. Non-redundant transcript assemblies generated after CD -176 HIT-EST clustering were translated into predicted coding sequences using TransDecoder v5.7.1 (Haas 177 et al. , 2013). The OrthoFinder v2.5 (Emms & Kelly, 2019) was used to infer orthogroups and 178 reconstruct a morph- level species tree. Maximum -likelihood gene trees were estimated for each 179 orthogroup and combined using the STAG algorithm to produce a genome -wide consensus topology, 180 which was rooted using STRIDE (Emms & Kelly, 2019). 181 Immediately after constructing the filtered gene- level expression matrix, we evaluated whether the 182 RWh morph exhibits intermediate expression consistent with simple parental inheritance or a 183 transgressive regulatory state relative to the two morphs observed in sympatry (PP and RW) using a 184 hybrid index (HI). Gene -level expression values were matched to morph identity, and genes with low 185 expression (mean ≤ 1 across samples) were removed. For each gene, g, mean expression was 186 calculated for each morph, and a hybrid index was computed along the parental PP–RW axis as: 187 𝐻𝐻𝐻𝐻𝐻𝐻= 𝜇𝜇𝜇𝜇𝜇𝜇ℎ, 𝐻𝐻− 𝜇𝜇𝜇𝜇𝜇𝜇, 𝐻𝐻 𝜇𝜇𝜇𝜇 𝜇𝜇, 𝐻𝐻− 𝜇𝜇𝜇𝜇𝜇𝜇, 𝐻𝐻 188 Where, 𝜇𝜇𝑚𝑚,𝐻𝐻 is the mean expression of gene g in morph m. This formula follows hybrid index 189 approaches used to quantify genomic or expression inheritance along parental gradients (Landry et al., 190 2007; McManus et al. , 2010). Genes with 𝜇𝜇RW,𝐻𝐻 = 𝜇𝜇PP,𝐻𝐻 were excluded for that gene because the 191 denominator was zero. Values between 0 and 1 indicate expression within the parental range, whereas 192 values 1 indicate transgressive expression outside the parental range. To account for differences 193 at the gene expression scale, the analysis was repeated after gene- wise standardization (z-scoring), 194 allowing detection of transgressiveness in relative expression profiles rather than absolute abundance. 195 Genes were classified into inheritance categories based on hybrid index thresholds: PP -like (0 ≤ HI ≤ 196 0.25), intermediate (0.25 < HI < 0.75), RW-like (0.75 ≤ HI ≤ 1), transgressive below PP ( HI 1). Genome -wide distributions of hybrid index values were visualized 198 using histograms, and category proportions were calculated to quantify the relative contribution of 199 additive, dominant, and transgressive regulatory patterns. 200 T ranscriptomic analyses and regulatory architecture 201 Following the early assessment of genome -wide expression inheritance and transgressiveness, we 202 conducted comprehensive analyses to characterize morph -associated regulatory divergence and the 203 transcriptional architecture of the RWh lineage. Gene- level differential expression among floral 204 morphs were assessed using the R -package DESeq2 v1.42 (Love et al. , 2014), with morph identity 205 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint specified as the focal factor and population included as a blocking factor, where appropriate. Count 206 data were modeled using a negative binomial framework, and significance was evaluated using Wald 207 tests with Benjamini–Hochberg false discovery rate (FDR) correction. Differentially expressed genes 208 provided an initial measure of morph- associated regulatory divergence and informed downstream 209 functional and network analyses. 210 To identify coordinated regulatory programs, weighted gene co -expression networks were 211 constructed using the R -package WGCNA v1.72-5 (Langfelder & Horvath, 2008) using TMM -212 normalized, log2-transformed floral expression data. An unsigned network was inferred using Pearson 213 correlations, with soft -thresholding parameters chosen to approximate a scale -free topology. Co-214 expression modules were identified based on topological overlap and defined using dynamic tree 215 cutting. Module eigengenes (first principal component of module expression) were used as low -216 dimensional summaries of regulatory activity. Module eigengenes formed the basis for all downstream 217 eco-evolutionary analyses. 218 Functional interpretation of transcriptional divergence was performed at the gene level using the R-219 package clusterProfiler v4.0 (Yu et al. , 2012), testing Gene Ontology (GO) and KEGG pathway 220 enrichment (Kanehisa, 2000) for differentially expressed genes and WGCNA modules. Enrichment 221 significance was assessed using Benjamini –Hochberg FDR correction. Analyses focused on pathways 222 associated with floral pigmentation, secondary metabolism, plastid function, and regulatory processes 223 to facilitate biological interpretation of morph-associated and RWh-specific regulatory patterns. 224 Partitioning morph-associated divergence versus geographic structure 225 To evaluate whether transcriptomic variation aligns more strongly with floral phenotype than with 226 spatial structure, module eigengenes were averaged within populations to generate population- level 227 regulatory profiles. Euclidean distance matrices were analyzed using permutational multivariate 228 analysis of variance (PERMANOVA) implemented in the adonis2 function in the R -package vegan 229 v2.6-4 (Oksanen et al., 2019) using 999 permutations, with morph identity as the explanatory factor. 230 Spatial structure was modeled using spatial eigenvectors derived from principal coordinates of 231 neighbor matrices (PCNM) calculated from geographic distances among populations (Borcard et al., 232 2018). 233 Distance-based redundancy analysis (dbRDA) was then performed using the capscale function in 234 the R-package vegan v2.6-4 (Oksanen et al., 2019), with morph identity as the explanatory factor and 235 spatial eigenvectors included as conditional variables. Morph identity was included as the focal 236 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint predictor, with spatial eigenvectors specified as conditional variables. Significance of constrained 237 models and marginal effects was assessed using permutation tests. 238 Variance partitioning analysis was then used to quantify the unique and shared contributions of 239 morph identity and geographic structure to transcriptomic divergence (Borcard et al., 1992; Peres-Neto 240 et al. , 2006). Adjusted R2 values were reported to avoid inflation of explained variance. A greater 241 proportion of variance uniquely attributable to morph identity relative to geography was interpreted as 242 evidence that transcriptional divergence is structured primarily by phenotype -associated selection 243 rather than neutral isolation by distance alone, consistent with recent eco- evolutionary genomic studies 244 using multivariate expression data (Forester et al., 2018; Capblancq et al., 2020). 245 Quantifying transcriptional canalization using within-population variance 246 Transcriptional canalization was quantified as reduced within -population variance of regulatory 247 modules (Waddington, 1942; Gibson & Wagner, 2000). For each population and morph, within-248 population variance was calculated as the standard deviation (SD) of module eigengene expression 249 across biological replicates. Module -level SDs were calculated separately for each population and 250 morph using workflows implemented in the R -package tidyverse v2.0.0 (Wickham et al. , 2019) . 251 Populations with only one replicate (n=1 in population S5; Table S1) were excluded from variance 252 estimation but retained for population- mean analyses. Low within- population variance for a given 253 module was interpreted as stronger regulatory constraint, consistent with stabilizing selection acting on 254 gene regulatory networks (Rifkin et al., 2005; Landry et al., 2007; Masel & Siegal, 2009). 255 Distributions of module-level variance were compared among floral morphs using Wilcoxon rank-256 sum tests, with false discovery rate (FDR) correction applied using the Benjamini–Hochberg procedure 257 (Benjamini & Hochberg, 1995). Reduced variance was interpreted as stronger regulatory constraint 258 consistent with stabilizing selection (Landry et al., 2007; McManus et al. , 2010). To quantify hybrid 259 destabilization, we calculated a hybrid variance inflation index for each module as: 260 𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑉= 𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑉ℎ 𝑚𝑚𝑚𝑚𝑉𝑉 𝑚𝑚 (𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑉, 𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑉) 261 Values greater than one indicate reduced canalization relative to parental morphs, and potential 262 regulatory incompatibility (Rieseberg et al. , 1999; Mack & Nachman, 2017). Modules exhibiting 263 elevated hybrid variance inflation were treated as candidates for regulatory destabilization and 264 potential phenotypic novelty. 265 Classifying hybrid inheritance at the regulatory module level 266 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint Regulatory inheritance was evaluated at the level of co -expression modules using population- mean 267 module eigengenes for PP, RW, and RWh. Inheritance categories were defined as additive, dominant 268 (PP-like or RW-like), or transgressive following established frameworks for RWh expression analysis 269 (Landry et al., 2007; McManus et al. , 2010). To minimize arbitrary thresholds, classification tolerance 270 was defined using module-specific parental variance, calculated as half of the pooled parental standard 271 deviation. This variance -informed approach aligns inheritance inference with the natural scale of 272 regulatory variation for each module and reduces misclassification due to stochastic noise (Romero et 273 al., 2012) . Inheritance categories were summarized across modules and compared with variance 274 inflation indices using Wilcoxon tests. 275

276 Discrete floral morphs show stable geographic occurrence and functional differentiation 277 Field sampling across Southwest China identified four recurrent floral morphs of Stellera 278 chamaejasme, pure pink (PP), red–white (RW), yellow–white (YW), pure yellow (YY), and a naturally 279 occurring red–white mosaic morph (RWh) occurring in zones of sympatry between PP and RW 280 populations (Fig. 2A; Supplementary Table S1). Each parental morph formed morphologically 281 consistent populations across geographically separated sites spanning broad elevational gradients 282 (2,676–4,210 m a.s.l.), indicating stable phenotypic recurrence across environmental and spatial 283 contexts (Fig. 2A). 284 Floral reflectance measurements revealed clear spectral differentiation among parental morphs (PP, 285 RW, YW, and YY) across the visible range (300–700 nm) (Fig. 2B). The PP morph exhibited low 286 reflectance at shorter wavelengths with a pronounced increase beyond ~600 nm, consistent with 287 anthocyanin-dominated pigmentation. In contrast, YW and YY morphs showed elevated reflectance 288 between ~500–600 nm and high reflectance at longer wavelengths, consistent with reduced 289 anthocyanin and increased carotenoid contributions . The RW morph displayed intermediate spectral 290 characteristics, with elevated reflectance in the green –yellow region. These distinct and non-291 overlapping spectral profiles support classification of morphs as discrete phenotypic categories rather 292 than continuous variation. The RWh morph was rare and restricted to sympatric zones between PP and 293 RW populations. Although spectral measurements were not obtained due to limited sampling, 294 representative inflorescences indicate a mosaic phenotype combining red floral tubes with mixed 295 white–pink lobes (Fig. 2C–G), distinguishing it from both parental morphs. 296 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint To ensure that downstream expression analyses were not confounded by assembly or annotation 297 bias, we evaluated transcriptome completeness and structural consistency across morphs and tissues. 298 BUSCO analyses indicated uniformly high completeness (97–99%) with comparable duplication and 299 fragmentation levels (Supplementary Fig. S1A). Expression- weighted contiguity (ExN50) increased 300 with transcript abundance and reached ~2.5 kb at Ex90 (Supplementary Fig. S1B), indicating high 301 assembly quality. Orthogroup analysis revealed extensive transcript sharing among morphs, indicating 302 a largely conserved gene repertoire ( Supplementary Fig. S1C). Mean–variance relationships followed 303 the expected negative binomial dispersion ( Supplementary Fig. S1D), and transcript length and 304 isoform distributions were comparable across datasets ( Supplementary Fig. S1E –K). RNA libraries 305 were of high quality (mapping rates ~95–97%, Q30 >94%), and Trinity assemblies showed consistent 306 contiguity and gene representation before and after filtering (Tables S2–S4). Together, these results 307 confirm that transcriptomic datasets are comparable across morphs and suitable for downstream 308 comparative analyses. 309 Phylogenetic placement and genome-wide expression inheritance of the RWh morph 310 To evaluate whether the mosaic phenotype of the RWh morph reflects intermediate evolutionary 311 ancestry, we reconstructed a morph- level phylogeny using genome -wide single -copy orthologs 312 inferred with OrthoFinder. The resulting topology was well resolved and did not support a simple 313 hybrid origin of RWh between the sympatric parental morphs (PP and RW) (Fig. 2H). Instead, RWh 314 was nested within the PP lineage as its closest relative, with RW forming a more distant sister clade. 315 This placement indicates that the mosaic phenotype of RWh does not correspond to phylogenetic 316 intermediary between PP and RW. Consistent with this phylogenetic pattern, orthogroup analysis 317 revealed that the vast majority of RWh genes were assigned to shared orthogroups (94%), with a low 318 proportion of morph- specific genes (3.7%) relative to parental morphs (Supplementary Table S5). 319 These results indicate that RWh does not represent a genomically distinct lineage but instead shares a 320 largely conserved gene repertoire with parental morphs. 321 We next assessed whether gene expression in RWh follows additive expectations along the PP –322 RW axis using a hybrid index framework. Genome -wide classification revealed that intermediate 323 expression was rare, representing only 6.89% of genes after scaling (Table 1). Parent -like expression 324 patterns were also uncommon (PP -like: 4.93%; RW -like: 2.16%). In contrast, the majority of genes 325 (86.02%) exhibited transgressive expression outside the parental range. This transgressive pattern was 326 asymmetric, with most genes showing expression below PP levels (74.33%), whereas a smaller 327 fraction exceeded RW expression (11.69%). 328 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint Together, the phylogenetic placement and genome -wide expression analyses indicate that RWh 329 does not represent a simple additive hybrid state. Instead, the mosaic phenotype is associated with 330 extensive transgressive regulatory variation and phylogenetic affinity to the PP lineage. These results 331 support a model in which RWh reflects a derived regulatory state characterized by widespread non-332 additive expression rather than intermediate inheritance. 333 Gene expression divergence reveals canalized regulatory states and transgressive expansion 334 Genome-wide gene expression profiles were strongly structured by floral morph rather than geographic 335 origin. Principal component analysis of normalized expression data separated individuals into discrete 336 clusters corresponding to the five morphs, with PC1 and PC2 explaining 21% and 11% of total 337 variance, respectively (Fig. 3A). Parental morphs formed compact, non- overlapping clusters, 338 indicating low within -morph transcriptional variance and consistent regulatory states across 339 populations. In contrast, the RWh morph occupied a displaced and expanded region of expression 340 space relative to the parental axis, consistent with increased regulatory dispersion rather than simple 341 additive inheritance. 342 Differential expression analyses revealed extensive morph -associated transcriptional divergence 343 (Fig. 3B–C; Supplementary Fig. S2). Across pairwise contrasts, thousands of genes were significantly 344 differentially expressed (FDR -adjusted P < 0.05), with both shared and morph- specific components , 345 with directionality dependent on the comparison order (i.e., relative up- or down-regulation between 346 morphs; Supplementary Fig. S2). Although a small core set of DEGs was shared among morphs, 347 substantial fractions were morph- specific, indicating that floral morphs correspond to discrete 348 transcriptional programs rather than continuous variation. 349 The PP morph showed the largest fraction of uniquely expressed genes (~12%), with consistent up-350 regulation of flavonoid and anthocyanin biosynthesis genes across multiple pairwise comparisons (e.g., 351 PP vs RW, PP vs RWh, PP vs YW; Supplementary Fig. S2). Key structural genes, including F3H, 352 ANS, CCR1, HCT, 4CL2, and L1AT , were repeatedly more highly expressed in PP relative to all other 353 morphs (Fig. 3C; Supplementary Fig. S2; Supplementary Table S6 ), consistent with elevated 354 anthocyanin flux underlying pink pigmentation. The PP morph also showed enrichment of regulatory 355 and stress -associated genes such as WRKY6, HST, PER12, and PRDX6 , suggesting coordinated 356 regulation of pigment biosynthesis and redox-related processes. The RW morph exhibited a smaller but 357 distinct set of DEGs (~7%), including BGLU18, CCD8, GT5, and HRPN, which are associated with 358 carotenoid cleavage, glycosylation, and secondary metabolite modification (Fig. 3B ; Supplementary 359 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint Table S6). Several phenylpropanoid-associated genes ( e.g., COMT, CTOMT1, BGLU24) were shared 360 between PP and RW, indicating partial overlap in regulatory architecture , as also reflected in their 361 pairwise contrasts (Supplementary Fig. S2). Yellow -based morphs (YW and YY) displayed strongly 362 divergent transcriptional profiles relative to PP and RW, with minimal overlap in DEGs across 363 pairwise comparisons (Supplementary Fig. S2). The YY morph uniquely expressed genes involved in 364 carotenoid metabolism and plastid function, including BGLU5, BGLU6, CCOMT, FOX2 , and PER47, 365 whereas the YW morph was characterized by genes such as BBX22, Z -ISO, F3GGT1 , and 366 CCOAOMT2 (Fig. 3B ; Supplementary Table S6 ). These patterns are consistent with pigmentation 367 primarily driven by carotenoid accumulation and reduced anthocyanin contribution. 368 The RWh morph exhibited a mixed but non- additive expression profile. While subsets of genes 369 displayed a distinct set of uniquely expressed genes (~11%), including BETA- OHASE, BGLU31, FGT, 370 SKIP1, PER10 , and TIFY6B , a substantial fraction displayed transgressive patterns relative to both 371 parental morphs (Fig. 3B ; Supplementary Table S6). These genes are associated with hormone 372 signaling, transcriptional regulation, and stress responses, indicating activation of regulatory pathways 373 not predominant in either parent. Notably, several core anthocyanin biosynthesis genes (e.g., ANS, 374 F3H) showed intermediate expression, whereas others (e.g., COMT, CTOMT1 ) were reduced relative 375 to both parents. In contrast, genes associated with stress response and signaling, including MSRB1 , 376 GC1, and NCED1 , were frequently up -regulated beyond parental levels, indicating activation of 377 regulatory pathways not predominant in either parental morph. 378 Hierarchical clustering of pigment - and metabolism -associated genes further highlighted 379 differences in regulatory stability among morphs ( Supplementary Fig. S3). Parental morphs displayed 380 highly coordinated expression patterns. Specifically, PP and RW consistently showed upregulation of 381 key phenylpropanoid and flavonoid pathway genes (e.g., HCT, BGLU13, L1AT, HST ), whereas YW 382 and YY exhibited coordinated downregulation of these genes, accompanied by increased expression of 383 carotenoid-associated genes such as Z- ISO and F3GGT1. In contrast, the RWh morph showed 384 heterogeneous expression across these clusters, combining intermediate ( ANS, F3H ), parent -biased, 385 and transgressive patterns (e.g., COMT, PER12, BGLU24 ). These results indicate that while parental 386 morphs maintain coordinated regulatory programs, the RWh morph reflects partial disruption of these 387 co-regulated gene networks. 388 Functional enrichment reflects metabolic bases of morph differentiation 389 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint Gene Ontology (GO) enrichment analyses revealed functional categories consistent with the metabolic 390 and regulatory basis of floral pigmentation (Supplementary Fig. S4). Across morphs, enriched 391 biological processes included phenylpropanoid and secondary metabolite biosynthesis, oxidation–392 reduction processes, and responses to stimuli, reflecting the coordinated metabolic and regulatory 393 demands of pigment production. Enrichment of plastid- and membrane-associated components further 394 highlights the cellular context of carotenoid biosynthesis and metabolite transport. 395 Morph-specific enrichment patterns were consistent with distinct pigmentation strategies. The PP 396 morph showed strong enrichment of flavonoid and phenylpropanoid pathways together with redox-397 related processes (Supplementary Fig. S4C), consistent with high anthocyanin flux and associated 398 oxidative regulation. The RW morph was enriched for regulatory and signaling categories 399 (Supplementary Fig. S4B), suggesting modulation of pigment pathways rather than full pathway 400 activation. In contrast, yellow -based morphs (YW and YY) were enriched for plastid -associated 401 metabolic processes (Supplementary Fig. S4E, F), consistent with carotenoid biosynthesis and storage 402 within plastids. The RWh morph exhibited broader and less specialized enrichment across metabolic, 403 regulatory, and stress-response categories (Supplementary Fig. S4D), indicating a more heterogeneous 404 transcriptional state. This contrasts with the more functionally coherent enrichment profiles observed 405 in parental morphs and is consistent with partial disruption of coordinated regulatory programs. 406 The KEGG pathway analyses further supported these patterns by highlighting differences in central 407 metabolic allocation (Fig. 3D; Supplementary Fig. S5). Core pathways such as carbon metabolism, 408 amino acid biosynthesis, and MAPK signaling were enriched across morphs, reflecting the metabolic 409 cost and regulatory control of pigment production. However, pathway representation differed among 410 morphs: PP and RW showed stronger enrichment for secondary metabolism and signaling pathways, 411 whereas YW and YY were enriched for lipid metabolism and plastid -associated processes. Notably, 412 the RWh morph showed the strongest enrichment of central metabolic pathways, including 413 glycolysis/gluconeogenesis and fatty acid metabolism, suggesting reorganization of primary metabolic 414 activity. Together, these results indicate that floral color polymorphism is embedded within broader 415 metabolic and regulatory networks rather than driven by isolated pigment pathways. 416 Modular regulatory architecture underlies floral morph differentiation 417 Weighted gene co -expression network analysis revealed a modular regulatory architecture underlying 418 floral morph differentiation (Fig. 4A; Supplementary Fig. S6). Genes clustered into discrete co -419 expression modules of varying sizes, each representing coordinated transcriptional programs. The 420 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint network construction achieved approximate scale -free topology at the selected soft -thresholding 421 power, with a clear trade -off between scale-free model fit and mean connectivity ( Supplementary Fig. 422 S6A, B). Hierarchical clustering of genes based on topological overlap delineated well -defined 423 modules of varying sizes, each summarized by a module eigengene ( Supplementary Fig. S6C). These 424 modules provide a systems -level framework for evaluating how gene expression is organized across 425 morphs. 426 Module–trait correlation analyses identified strong and morph -specific regulatory associations 427 (Supplementary Fig. S6D). Several modules were positively correlated with individual parental 428 morphs, indicating that each morph is characterized by a distinct combination of regulatory programs 429 rather than uniform shifts in gene expression. In contrast, modules associated with the RWh morph did 430 not simply mirror parental patterns, suggesting that its regulatory organization is not intermediate at 431 the network level. Hierarchical clustering of module eigengenes further revealed higher -order structure 432 among modules, grouping regulatory programs with shared expression dynamics (Supplementary Fig. 433 S6E). Sample clustering based on module eigengenes confirmed clear separation among morphs 434 without evidence of outliers (Supplementary Fig. S6F), supporting the robustness of the inferred 435 network structure. Functional enrichment of morph- associated modules indicated that regulatory 436 architecture reflects biological specialization. Modules associated with PP and RW were enriched for 437 phenylpropanoid and flavonoid biosynthesis pathways, whereas modules associated with YW and YY 438 were enriched for carotenoid metabolism and plastid -associated processes (Supplementary Fig. S7). 439 Modules associated with the RWh morph showed enrichment for regulatory and stress -response 440 functions, consistent with altered transcriptional states relative to parental morphs. The absence of 441 intermediate module -level patterns in RWh further indicates that regulatory variation in this morph 442 reflects reorganization of network structure rather than additive inheritance. 443 Canalization and selective hybrid destabilization shape regulatory divergence 444 To disentangle morph- associated regulatory divergence from neutral spatial structure, we partitioned 445 variation in population- mean module eigengenes using PERMANOVA, distance -based redundancy 446 analysis (dbRDA), and variance partitioning (Table 2). Canonical analysis of principal coordinates 447 based on population-mean module eigengenes further separated morphs along major axes of regulatory 448 variation, with the RWh morph positioned between PP and RW but shifted away from the direct 449 parental axis, occupying a distinct region of expression space rather than a simple mid- parent position 450 (Fig. 4B). Morph identity explained a substantial fraction of transcriptomic variation (R² = 0.45, P = 451 0.001), indicating strong phenotype -associated structuring of regulatory profiles. After accounting for 452 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint spatial structure using geographic eigenvectors, morph identity retained a marginal but persistent effect 453 (dbRDA: F = 1.43, P = 0.08), suggesting that morph- associated divergence is not fully explained by 454 geographic distance alone (Fig. 4B; Table 2). Variance partitioning further revealed a larger unique 455 contribution of morph identity (adjusted R² = 0.16) relative to geography (adjusted R² = 0.12), with a 456 smaller shared component (adjusted R² = 0.05). Together, these results indicate that transcriptomic 457 divergence aligns more strongly with morph identity than with spatial structure (Fig. 4C). 458 Within-population variance of module eigengenes provided a quantitative measure of 459 transcriptional constraint across morphs (Fig. 4D; Supplementary Fig. S8A). Parental morphs (PP, 460 RW, YW, YY) consistently exhibited low variance across modules, indicating stable and tightly 461 regulated expression states. In contrast, the RWh morph showed a similar median variance but a 462 pronounced tail of high -variance modules, indicating elevated variability in a subset of regulatory 463 programs. 464 Hybrid variance inflation analysis further quantified this pattern. The median inflation index across 465 modules was low (0.57), indicating that most regulatory modules remained stable in RWh relative to 466 parental morphs (Table 2). However, a small subset of modules exhibited extreme inflation, with 467 values exceeding 25 in modules such as MEskyblue3 and MEmagenta (Table 3). These modules 468 (arbitrarily labeled by WGCNA) represent distinct gene regulatory subnetwork, suggesting that 469 regulatory instability is localized to specific functional pathways rather than diffusely distributed 470 across the transcriptome (Fig. 4E; Table 2). Module- level inheritance analysis revealed structured but 471 heterogeneous regulatory outcomes in RWh (Fig. 4F; Supplementary Fig. S8B; Table 3). Across all 472 modules (n = 53), 26.4% (14 modules) exhibited additive inheritance, while 35.8% (28.3% dominant -473 PP and 7.5% dominant-RW) showed dominant patterns biased toward one parental morph. In contrast, 474 37.7% of modules (20 modules) displayed transgressive inheritance, with expression extending beyond 475 the parental range, consistent with the emergence of novel regulatory states (Fig. 4F; Table 2). 476 Modules exhibiting transgressive inheritance tended to show higher variance inflation than additive 477 or dominant modules, although this relationship was not statistically significant (Wilcoxon P = 0.25; 478 Table 2). The most extreme cases of hybrid regulatory disruption were confined to a few modules; for 479 example, MEskyblue3 showed a hybrid variance inflation index of 31.52 together with transgressive 480 inheritance. Together, these results indicate that most gene regulatory modules remain stable across 481 morphs, whereas a subset is selectively destabilized in the RWh morph. This pattern supports a model 482 in which transcriptional systems are largely robust but permit localized breakdown of regulatory 483 constraint, generating transgressive expression and potential phenotypic novelty. 484 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint

485 Morph identity defines stable transcriptional states across geographic space 486 The central goal of this study was to determine whether flower -color polymorphism in Stellera 487 chamaejasme reflects transient phenotypic variation or stable regulatory states associated with morph 488 identity. Conceptually, transient phenotypic variation would be expected to produce weak and 489 inconsistent transcriptional differentiation among morphs, with expression patterns primarily 490 structured within populations or across environmental gradients. In contrast, persistent morph-491 associated regulatory states would be expected to show stronger and more consistent transcriptional 492 divergence among morphs than within morphs or populations, remain stable across geographic space, 493 and persist after accounting for spatial structure. Across populations spanning broad elevational 494 gradients (2,676–4,210 m a.s.l.) in Southwest China, transcriptomic variation was consistently 495 structured by morph rather than geographic distance. Partitioning analyses further showed that morph 496 identity explained a larger and distinct component of expression variation than geography, and this 497 signal remained after explicit spatial correction. These results indicate that floral morphs correspond to 498 reproducible phenotype-associated transcriptional programs and reveal an isolation-by- morph pattern 499 at the regulatory level. 500 Disentangling phenotype-associated divergence from spatial structure is particularly challenging in 501 topographically complex regions such as the Qinghai –Tibetan Plateau, where environmental and 502 geographic gradients are strongly correlated (Meirmans, 2012; Wang & Bradburd, 2014). By explicitly 503 modeling spatial effects, our analyses provide a conservative test of morph- associated divergence and 504 reduce the likelihood of spurious inference. The persistence of morph- associated transcriptional 505 structure after spatial correction indicates that regulatory differentiation is unlikely to arise solely from 506 isolation by distance or demographic history. 507 Although transcriptomic data alone cannot directly demonstrate the mode of selection, the 508 observed patterns are consistent with morph- specific regulatory configurations potentially maintained 509 by phenotype -associated or selective processes. Comparable phenotype -linked regulatory 510 differentiation has been reported in several floral polymorphism systems, including monkeyflowers 511 Mimulus L. (Schemske & Bradshaw, 1999)/Mimulus sect. Erythranthe (Spach) Greene (Wenzell et al., 512 2025), Penstemon Schmidel (Wessinger et al. , 2014) , and Phlox L. (Hopkins & Rausher, 2012) , as 513 well as in eco -genomic studies showing that phenotypic predictors can explain expression or genomic 514 variation beyond spatial structure (Wang & Bradburd, 2014; Forester et al. , 2016; Capblancq et al. , 515 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint 2020). The shared morph–geography component observed here further suggests that environmental 516 gradients may also influence morph distributions through variation in abiotic conditions, pollinator 517 communities, or resource availability. Such coupling is expected when selection varies along 518 environmental gradients but favors discrete alternative phenotypes rather than continuous variation 519 (Schemske & Bradshaw, 1999; Hopkins & Rausher, 2012). Importantly, the substantial morph-specific 520 component independent of geography suggests that regulatory divergence is associated with morph 521 identity itself rather than arising solely as a by -product of environmental variation. Collectively, these 522 findings indicate that flower -color morphs in S. chamaejasme represent stable regions of 523 transcriptional state space that recur across heterogeneous landscapes. 524 Canalization of regulatory networks underlies morph stability 525 Canalization is a key feature of stabilizing selection, where phenotypic and gene expression variation 526 is reduced around an optimal regulatory state (Waddington, 1942; Gibson & Wagner, 2000; Masel & 527 Siegal, 2009). Consistent with this expectation, parental morphs exhibited low within- population and 528 within-morph variance in co -expression module eigengenes, even among geographically separated 529 populations spanning approximately 23–410 km . Notably, this pattern is observed at the level of 530 regulatory modules rather than individual genes, suggesting that morph differentiation is embedded 531 within network-level expression architectures. Such module -level stability reflects regulatory systems 532 that buffer expression against environmental and genetic perturbation. 533 Module-level canalization reflects regulatory systems that buffer expression against environmental 534 and genetic perturbation (Rifkin et al. , 2005; Hodgins -Davis & Townsend, 2009; McManus et al. , 535 2010). The persistence of these low -variance expression regimes across geographically separated 536 populations suggests that each floral morph occupies a distinct and stable region of transcriptional 537 space consistent with maintenance by selective or ecological filtering processes . Independent 538 population genomic evidence showing strong regional structure and limited gene flow further indicates 539 that these regulatory configurations persist despite demographic differentiation (Rana et al., 2024). 540 This pattern was particularly evident in pigmentation -associated regulatory modules . Genes 541 involved in phenylpropanoid, flavonoid, and carotenoid metabolism were consistently organized into 542 morph-specific co -expression modules exhibiting low within-morph variance. Because floral 543 coloration depends on coordinated flux through multi -step biosynthetic pathways, robustness is 544 expected to emerge at the pathway and network levels rather than at individual loci (Grotewold, 2006; 545 Davies et al., 2012; Albert et al. , 2014). The stability of these pigment-associated regulatory modules 546 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint across populations therefore provides a direct mechanistic link between transcriptional canalization 547 and the persistence of discrete flower-color phenotypes, reinforcing the view that morph differentiation 548 in S. chamaejasme is deeply embedded in regulatory architecture rather than driven by isolated 549 pigment genes (Rausher, 2008; Yuan et al., 2013; Monniaux, 2023). 550 Transgressive morphs reveal limits of regulatory stability 551 The naturally occurring RWh morph provides insight into the limits of these regulatory states. 552 Although field observations suggested a mosaic phenotype in zones of sympatry between PP and RW 553 populations, ortholog- based phylogenetic and expression analyses do not support a simple additive 554 hybrid interpretation. Instead, genome -wide hybrid index analyses revealed pervasive non -additive 555 regulation, with approximately 86% of genes exhibiting expression outside the parental range. These 556

indicate that RWh represents a transgressive regulatory morphotype rather than an intermediate 557 expression state. 558 When parental morphs correspond to stabilized regulatory configurations, recombination or 559 regulatory reorganization may disrupt co-adapted interactions among network components, exposing 560 epistatic interactions and relaxing stabilizing constraints (Muller, 1942; Dobzhansky, 1982; Rieseberg 561 et al., 1999). Consistent with this pattern, the RWh morph exhibited elevated expression variance in a 562 subset of co- expression modules. Importantly, this destabilization was modular rather than genome -563 wide. Most modules retained parental -like stability and predictable inheritance, whereas a limited 564 subset showed pronounced variance inflation and transgressive expression. This selective 565 destabilization aligns with theoretical and empirical evidence that regulatory robustness is unevenly 566 distributed across gene networks, where core modules remain buffered while others are more sensitive 567 to perturbation (Gibson & Wagner, 2000; Wagner, 2008; Paaby & Rockman, 2014), and is further 568 supported by pathway-specific hybrid misexpression across diverse taxa (McManus et al., 2010; Mack 569 & Nachman, 2017; Runemark et al., 2025). 570 Module-level inheritance patterns further indicate that regulatory novelty arises through localized 571 disruption of network organization. While many modules exhibited additive or dominant expressions, a 572 substantial proportion showed transgressive regulation beyond the parental range. Such expression 573 patterns are commonly interpreted as signatures of epistatic interactions, compensatory regulatory 574 evolution, or breakdown of co- adapted regulatory complexes (Rieseberg et al. , 1999; Stelkens & 575 Seehausen, 2009; Camper et al. , 2025; Runemark et al. , 2025). In S. chamaejasme, the coexistence of 576 stable regulatory modules alongside a subset of transgressive modules suggests that gene regulatory 577 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint networks can remain globally robust while permitting localized release of expression variation. 578 Notably, several of the modules exhibiting transgressive regulation were enriched for pigment 579 biosynthesis and secondary metabolism, directly linking reduced regulatory constraint to floral color 580 variation. 581 Regulatory integration links flower pigmentation to broader physiological organization 582 Our results further indicate that flower -color morphs in S. chamaejasme are not defined solely by 583 pigment biosynthesis but correspond to coordinated transcriptional states involving multiple 584 physiological processes. Functional enrichment analyses revealed that morph- associated expression 585 divergence extends beyond anthocyanin and flavonoid pathways to include carbon metabolism, redox 586 processes, hormone signaling, and developmental regulation. These findings indicate that flower -color 587 polymorphisms reflect integrated physiological states rather than isolated biochemical traits. 588 This integrative organization is increasingly recognized as a defining feature of floral trait 589 evolution. Pigment pathways are tightly linked to primary metabolism, stress responses, and 590 developmental regulation (Smith & Rausher, 2011; Martínez -Harms et al. , 2022; Yeo & Moyroud, 591 2025). Anthocyanins, for example, contribute not only to pollinator attraction but also to UV 592 protection, thermal regulation, and oxidative stress mitigation, particularly in high- elevation 593 environments (Neill & Gould, 2003; Davies et al. , 2012; Landi et al. , 2015). Because these functions 594 share metabolic substrates and regulatory control, selection acting on ecological performance may 595 indirectly contribute to stabilizing pigment-associated regulatory programs (Paaby & Rockman, 2014; 596 Wessinger & Rausher, 2015). In this context, module-level canalization may contribute to maintaining 597 functional integration by ensuring that pigmentation remains coordinated with broader physiological 598 processes (Gibson & Wagner, 2000; Masel & Siegal, 2009; Martínez -Harms et al. , 2022) , whereas 599 disruption of these modules can generate both expression variability and phenotypic novelty. The 600 transgressive RWh morph illustrates this pattern, showing altered regulation in modules enriched for 601 metabolic and regulatory functions, suggesting partial decoupling of coordinated gene networks 602 (Landry et al., 2007; McManus et al., 2010; Mack & Nachman, 2017). In S. chamaejasme, this pattern 603 implies that regulatory integration both stabilizes discrete floral phenotypes and constrains the range of 604 viable transgressive expression states. 605 Evolutionary implications and future directions 606 Together, our results indicate that flower -color morphs in S. chamaejasme correspond to discrete and 607 repeatable transcriptional states characterized by low within -morph variance and strong phenotype -608 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint associated structure across populations. These patterns are consistent with morph- specific regulatory 609 configurations potentially maintained by stabilizing or environmentally mediated processes rather than 610 neutral demographic structure alone. The presence of extensive transgressive expression and localized 611 variance inflation in the RWh morph further suggests that recombination or regulatory reorganization 612 may expose incompatibility among co- adapted network components, generating phenotypic novelty 613 through localized disruption of regulatory buffering (Seehausen et al., 2014). 614 Several limitations expands from considering sampling strategy to the fitness assessment . 615 Transcriptomic analyses were restricted to a specific tissue and developmental stage, and the strength 616 of canalization may vary across ontogeny or environmental conditions. Sampling of the transgressive 617 morph was limited due to its rarity, which may reduce power to detect subtle effects. In addition, while 618 the observed patterns are consistent with regulatory constraint, direct tests of selection and fitness 619 consequences will be required to establish the evolutionary mechanisms underlying morph persistence. 620 Future work linking regulatory architecture to ecological performance will be essential for 621 evaluating the evolutionary consequences of these transcriptional states. Experimental approaches such 622 as pollinator preference assays, reciprocal transplant experiments, and functional validation of 623 regulatory loci could clarify the fitness consequences of alternative morph- specific regulatory 624 programs. Integrating transcriptomic architecture with population genomic and ecological data will be 625 necessary to understand how regulatory robustness shapes the persistence, breakdown, and 626 evolutionary potential of polymorphic traits in natural populations. 627

628 We thank Zhen -Yu Lv for assistance with field sampling, Dr. Zhe Chen for taking the reflectance 629 measurements, and Ms. Song Min- shu for laboratory support. This work was supported by the Key 630 Projects of the Joint Fund of the National Natural Science Foundation of China (U23A20149), Yunnan 631 Key R&D Program (202403AC00028), the Second Tibetan Plateau Scientific Expedition and Research 632 (STEP) Program (2024QZKK0200), the Research Fund for International Young Scientists of the 633 National Natural Science Foundation of China (32150410356), the National Natural Science 634 Foundation of China (32322006), and the Young Academic and Technical Leader Raising Foundation 635 of Yunnan Province (2019HB039). The first author was additionally supported by the CAS President’s 636 International Fellowship Initiative (PIFI) postdoctoral fellowship (2021PB0034). 637 Competing interests 638 The authors declare no competing interests. 639 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint Author contributions 640 S.K.R., T.D. and H.S. planned and designed the research. S.K.R., H.K.R. and C.J. performed the 641 research and conducted field and laboratory work. S.K.R. and H.K.R. performed transcriptome 642 assembly and analysed the data. S.K.R. conducted eco- evolutionary and network analyses. S.K.R. and 643 H.K.R. prepared the figures. T.D. and H.S. supervised the project and acquired funding. S.K.R. wrote 644 the manuscript, and all authors contributed to reviewing, editing, and approving the final version. 645 Data availability 646 Raw RNA- seq reads generated in this study were deposited in the NCBI Sequence Read Archive 647 (SRA) under BioProject accession PRJNA1451939. Processed data, including expression matrices, 648 differential expression results, co-expression modules, and functional annotation outputs, are available 649 in the GitHub repository: https://github.com/santoshkumarrana/stellera -regulatory-network-650 polymorphism and archived in Zenodo ( 10.5281/zenodo.19307586) (Rana, 2026) . All data will be 651 publicly accessible upon publication. 652 Code availability 653 Custom scripts developed for this manuscript are publicly available on GitHub: 654 https://github.com/santoshkumarrana/Canalized-regulatory-networks-stabilize-floral-polymorphism-655 and-enable-localized-transgression. A version of the repository has been archived at Zenodo 656 (10.5281/zenodo.19307586) (Rana, 2026) to ensure long-term accessibility. 657 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint

658 Albert NW, Davies KM, Schwinn KE. 2014. Gene regulation networks generate diverse 659 pigmentation patterns in plants. Plant Signaling & Behavior 9: e29526. 660 Benjamini Y, Hochberg Y. 1995. Controlling the false discovery rate: A practical and powerful 661 approach to multiple testing. Journal of the Royal Statistical Society Series B: Statistical 662 Methodology 57: 289–300. 663 Bolger AM, Lohse M, Usadel B. 2014. Trimmomatic: a flexible trimmer for Illumina sequence 664 data. Bioinformatics 30: 2114–2120. 665 Borcard D, Gillet F, Legendre P. 2018. Numerical Ecology with R. Cham: Springer 666 International Publishing. 667 Borcard D, Legendre P, Drapeau P. 1992. Partialling out the spatial component of ecological 668 variation. Ecology 73: 1045–1055. 669 Boyle EA, Li YI, Pritchard JK. 2017. An expanded view of complex traits: From polygenic to 670 omnigenic. Cell 169: 1177–1186. 671 Camacho C, Coulouris G, Avagyan V, Ma N, Papadopoulos J, Bealer K, Madden TL. 2009. 672 BLAST+: Architecture and applications. BMC Bioinformatics 10: 421. 673 Camper BT, Kanes AS, Laughlin ZT, Manuel RT, Bewick SA. 2025. Transgressive hybrids 674 as hopeful holobionts. Microbiome 13: 19. 675 Cantalapiedra CP, Hernández-Plaza A, Letunic I, Bork P, Huerta-Cepas J. 2021. eggNOG-676 mapper v2: Functional annotation, orthology assignments, and domain prediction at the 677 metagenomic scale (K Tamura, Ed.). Molecular Biology and Evolution 38: 5825–5829. 678 Capblancq T, Fitzpatrick MC, Bay RA, Exposito-Alonso M, Keller SR. 2020. Genomic 679 prediction of (mal)adaptation across current and future climatic landscapes. Annual 680 Review of Ecology, Evolution, and Systematics 51: 245–269. 681 Chen H, Berg CS, Samuli M, Sotola VA, Sweigart AL, Yuan Y, Fishman L. 2025. The 682 genetic architecture of floral trait divergence between hummingbird‐ and self‐pollinated 683 monkeyflower (Mimulus) species. New Phytologist 245: 2255–2267. 684 Dauphin B, Rellstab C, Wüest RO, Karger DN, Holderegger R, Gugerli F, Manel S. 2023. 685 Re-thinking the environment in landscape genomics. Trends in Ecology & Evolution 38: 686 261–274. 687 Davies KM, Albert NW, Schwinn KE. 2012. From landing lights to mimicry: the molecular 688 regulation of flower colouration and mechanisms for pigmentation patterning. Functional 689 Plant Biology 39: 619–638. 690 Dobzhansky T. 1982. Genetics and the origin of species. New York: Columbia University 691 Press. 692 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint Emms DM, Kelly S. 2019. OrthoFinder: Phylogenetic orthology inference for comparative 693 genomics. Genome Biology 20: 238. 694 Félix M-A, Barkoulas M. 2015. Pervasive robustness in biological systems. Nature Reviews 695 Genetics 16: 483–496. 696 Forester BR, Jones MR, Joost S, Landguth EL, Lasky JR. 2016. Detecting spatial genetic 697 signatures of local adaptation in heterogeneous landscapes. Molecular Ecology 25: 104–698 120. 699 Forester BR, Lasky JR, Wagner HH, Urban DL. 2018. Comparing methods for detecting 700 multilocus adaptation with multivariate genotype–environment associations. Molecular 701 Ecology 27: 2215–2233. 702 Fu L, Niu B, Zhu Z, Wu S, Li W. 2012. CD-HIT: Accelerated for clustering the next-703 generation sequencing data. Bioinformatics 28: 3150–3152. 704 Gibson G, Wagner G. 2000. Canalization in evolutionary genetics: A stabilizing theory? 705 BioEssays 22: 372–380. 706 Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, Adiconis X, Fan L, 707 Raychowdhury R, Zeng Q, et al. 2011. Full-length transcriptome assembly from RNA-708 Seq data without a reference genome. Nature Biotechnology 29: 644–652. 709 Gray SM, McKinnon JS. 2007. Linking color polymorphism maintenance and speciation. 710 Trends in Ecology & Evolution 22: 71–79. 711 Grotewold E. 2006. The genetics and biochemistry of floral pigments. Annual Review of Plant 712 Biology 57: 761–780. 713 Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, Couger MB, 714 Eccles D, Li B, Lieber M, et al. 2013. De novo transcript sequence reconstruction from 715 RNA-seq using the Trinity platform for reference generation and analysis. Nature 716 Protocols 8: 1494–1512. 717 Hodgins-Davis A, Townsend JP. 2009. Evolving gene expression: from G to E to G×E. Trends 718 in Ecology & Evolution 24: 649–658. 719 Hopkins R, Rausher MD. 2012. Pollinator-mediated selection on flower color allele drives 720 reinforcement. Science 335: 1090–1092. 721 Horta-Lacueva QJ-B, Jónsson ZO, Thorholludottir DAV, Hallgrímsson B, Kapralova KH. 722 2023. Rapid and biased evolution of canalization during adaptive divergence revealed by 723 dominance in gene expression variability during Arctic charr early development. 724 Communications Biology 6: 897. 725 Kadelka C. 2026. Canalization as a stabilizing principle of gene regulatory networks: A discrete 726 dynamical systems perspective. NPJ Systems Biology and Applications 12: 30. 727 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint Kadelka C, Murrugarra D. 2024. Canalization reduces the nonlinearity of regulation in 728 biological networks. NPJ Systems Biology and Applications 10: 67. 729 Kanehisa M. 2000. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Research 730 28: 27–30. 731 Kitano J, Kagawa K, Tsuchimatsu T, Yamaguchi R, Yamamichi M. 2025. The genomics of 732 discrete polymorphisms maintained by disruptive selection. Trends in Ecology & 733 Evolution 40: 1023–1034. 734 Landi M, Tattini M, Gould KS. 2015. Multiple functional roles of anthocyanins in plant-735 environment interactions. Environmental and Experimental Botany 119: 4–17. 736 Landry CR, Hartl DL, Ranz JM. 2007. Genome clashes in hybrids: insights from gene 737 expression. Heredity 99: 483–493. 738 Langfelder P, Horvath S. 2008. WGCNA: An R package for weighted correlation network 739 analysis. BMC Bioinformatics 9: 559. 740 Li J, Wang S, Miehe G, Opgenoorth L, Yang H, Wu D, Mao K. 2025. Ecological selection as 741 drivers during early speciation: Insights from two allopatric cypress species in the 742 Himalaya. Molecular Ecology 34: e70072. 743 Love MI, Huber W, Anders S. 2014. Moderated estimation of fold change and dispersion for 744 RNA-seq data with DESeq2. Genome Biology 15: 550. 745 Mack KL, Nachman MW. 2017. Gene regulation and speciation. Trends in Genetics 33: 68–80. 746 Maia R, Gruson H, Endler JA, White TE. 2019. PAVO 2: New tools for the spectral and spatial 747 analysis of colour in R (RB O’Hara, Ed.). Methods in Ecology and Evolution 10: 1097–748 1107. 749 Martínez-Harms J, Guerrero PC, Martínez-Harms MJ, Poblete N, González K, Stavenga 750 DG, Vorobyev M. 2022. Mechanisms of flower coloring and eco-evolutionary 751 implications of massive blooming events in the Atacama Desert. Frontiers in Ecology 752 and Evolution 10: 957318. 753 Masel J, Siegal ML. 2009. Robustness: mechanisms and consequences. Trends in Genetics 25: 754 395–403. 755 McManus CJ, Coolon JD, Duff MO, Eipper-Mains J, Graveley BR, Wittkopp PJ. 2010. 756 Regulatory divergence in Drosophila revealed by mRNA-seq. Genome Research 20: 757 816–825. 758 Meirmans PG. 2012. The trouble with isolation by distance. Molecular Ecology 21: 2839–2846. 759 Monniaux M. 2023. Unusual suspects in flower color evolution. Science 379: 534–535. 760 Muller HJ. 1942. Isolating mechanisms, evolution, and temperature. 761 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint Neill SO, Gould KS. 2003. Anthocyanins in leaves: Light attenuators or antioxidants? 762 Functional Plant Biology 30: 865–873. 763 Oksanen J, Blanchet FG, Kindt R, Legendre P, Friendly M, McGlinn D, Minchin PR, 764 O’Hara RB, Simpson GL, Solymos P, et al. 2019. vegan: Community ecology package. 765 Paaby AB, Rockman MV. 2014. Cryptic genetic variation: evolution’s hidden substrate. Nature 766 Reviews Genetics 15: 247–258. 767 Patro R, Duggal G, Love MI, Irizarry RA, Kingsford C. 2017. Salmon provides fast and bias-768 aware quantification of transcript expression. Nature Methods 14: 417–419. 769 Payne JL, Wagner A. 2019. The causes of evolvability and their evolution. Nature Reviews 770 Genetics 20: 24–38. 771 Peres-Neto PR, Legendre P, Dray S, Borcard D. 2006. Variation partitioning of species data 772 matrices: Estimation and comparison of fractions. Ecology 87: 2614–2625. 773 Rana SK. 2026. Canalized gene regulatory networks stabilize floral polymorphism and enable 774 modular transgressive expression in Stellera chamaejasme. v1.0.0, Zenodo, 29 Mar. 775 2026, https://doi.org/10.5281/zenodo.19307586 776 Rana SK, Rana HK, Landis JB, Kuang T, Chen J, Wang H, Deng T, Davis CC, Sun H. 777 2024. Pleistocene glaciation advances the cryptic speciation of Stellera chamaejasme L. 778 in a major biodiversity hotspot. Journal of Integrative Plant Biology 66: 1192–1205. 779 Rausher MD. 2008. Evolutionary transitions in floral color. International Journal of Plant 780 Sciences 169: 7–21. 781 Ravinet M, Faria R, Butlin RK, Galindo J, Bierne N, Rafajlović M, Noor MAF, Mehlig B, 782 Westram AM. 2017. Interpreting the genomic landscape of speciation: a road map for 783 finding barriers to gene flow. Journal of Evolutionary Biology 30: 1450–1477. 784 Rellstab C, Gugerli F, Eckert AJ, Hancock AM, Holderegger R. 2015. A practical guide to 785 environmental association analysis in landscape genomics. Molecular Ecology 24: 4348–786 4370. 787 Rieseberg LH, Archer MA, Wayne RK. 1999. Transgressive segregation, adaptation and 788 speciation. Heredity 83: 363–372. 789 Rifkin SA, Houle D, Kim J, White KP. 2005. A mutation accumulation assay reveals a broad 790 capacity for rapid evolution of gene expression. Nature 438: 220–223. 791 Robinson MD, McCarthy DJ, Smyth GK. 2010. edgeR: A Bioconductor package for 792 differential expression analysis of digital gene expression data. Bioinformatics 26: 139–793 140. 794 Romero IG, Ruvinsky I, Gilad Y. 2012. Comparative studies of gene expression and the 795 evolution of gene regulation. Nature Reviews Genetics 13: 505–516. 796 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint Runemark A, Moore EC, Larson EL. 2025. Hybridization and gene expression: Beyond 797 differentially expressed genes. Molecular Ecology 34: e17303. 798 Schemske DW, Bradshaw HD. 1999. Pollinator preference and the evolution of floral traits in 799 monkeyflowers (Mimulus). Proceedings of the National Academy of Sciences 96: 11910–800 11915. 801 Seehausen O, Butlin RK, Keller I, Wagner CE, Boughman JW, Hohenlohe PA, Peichel CL, 802 Saetre G-P, Bank C, Brännström Å, et al. 2014. Genomics and the origin of species. 803 Nature Reviews Genetics 15: 176–192. 804 Simão FA, Waterhouse RM, Ioannidis P, Kriventseva EV, Zdobnov EM. 2015. BUSCO: 805 assessing genome assembly and annotation completeness with single-copy orthologs. 806 Bioinformatics 31: 3210–3212. 807 Smith SD. 2016. Pleiotropy and the evolution of floral integration. New Phytologist 209: 80–85. 808 Smith SD, Rausher MD. 2011. Gene loss and parallel evolution contribute to species difference 809 in flower color. Molecular Biology and Evolution 28: 2799–2810. 810 Stelkens R, Seehausen O. 2009. Genetic distance between species predicts novel trait 811 expression in their hybrids. Evolution 63: 884–897. 812 Sun S-F, Li Y-B, Leng J, Zhang D-H, Zhang W-L, Zhang B. 2024. Multiple insect pollination 813 contributes to differential phenotypic selection on floral traits in Stellera chamaejasme. 814 BMC Plant Biology 24: 1141. 815 Tigano A, Friesen VL. 2016. Genomics of local adaptation with gene flow. Molecular Ecology 816 25: 2144–2164. 817 Waddington CH. 1942. Canalization of development and the inheritance of acquired characters. 818 Nature 150: 563–565. 819 Wagner A. 2008. Robustness and evolvability: A paradox resolved. Proceedings of the Royal 820 Society B: Biological Sciences 275: 91–100. 821 Wagner A. 2011. The origins of evolutionary innovations: a theory of transformative change in 822 living systems. Oxford: Oxford University Press. 823 Wang IJ, Bradburd GS. 2014. Isolation by environment. Molecular Ecology 23: 5649–5662. 824 Wenzell KE, Neequaye M, Paajanen P, Hill L, Brett P, Byers KJRP. 2025. Within-species 825 floral evolution reveals convergence in adaptive walks during incipient pollinator shift. 826 Nature Communications 16: 2721. 827 Wessinger CA, Hileman LC, Rausher MD. 2014. Identification of major quantitative trait loci 828 underlying floral pollination syndrome divergence in Penstemon. Philosophical 829 Transactions of the Royal Society B: Biological Sciences 369: 20130349. 830 Wessinger CA, Rausher MD. 2015. Ecological transition predictably associated with gene 831 degeneration. Molecular Biology and Evolution 32: 347–354. 832 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint Wickham H, Averick M, Bryan J, Chang W, McGowan L, François R, Grolemund G, 833 Hayes A, Henry L, Hester J, et al. 2019. Welcome to the Tidyverse. Journal of Open 834 Source Software 4: 1686. 835 Yeaman S, Whitlock MC. 2011. The genetic architecture of adaptation under migration-836 selection balance: The genetic architecture of local adaptation. Evolution 65: 1897–1911. 837 Yeo MTS, Moyroud E. 2025. Under the rainbow: Novel insights on the mechanisms driving the 838 development and evolution of petal pigmentation. Current Opinion in Plant Biology 86: 839 102743. 840 Yu G, Wang L-G, Han Y, He Q-Y. 2012. clusterProfiler: An R package for comparing 841 biological themes among gene clusters. OMICS: A Journal of Integrative Biology 16: 842 284–287. 843 Yuan Y-W, Byers KJ, Bradshaw H. 2013. The genetic control of flower–pollinator specificity. 844 Current Opinion in Plant Biology 16: 422–428. 845 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint Supporting Information (brief legends) 846 Supplementary Figures 847 Supplementary Fig. S1 Transcriptome assembly quality, expression characteristics, and 848 transcript diversity. 849 Supplementary Fig. S2 Pairwise differential gene expression among flower-color morphs. 850 Supplementary Fig. S3 Gene expression heatmap for pigment - and metabolism -associated 851 genes. 852 Supplementary Fig. S4 Gene Ontology enrichment of differentially expressed genes across 853 flower-color morphs. 854 Supplementary Fig. S5 KEGG pathway enrichment of differentially expressed genes among 855 flower-color morphs. 856 Supplementary Fig. S6 Weighted gene co -expression network analysis (WGCNA) identifies 857 modules associated with flower-color morphs. 858 Supplementary Fig. S7 Pathway representation and functional network structure of pigment -859 associated genes. 860 Supplementary Fig. S8 Within -morph expression variance and inheritance patterns of co-861 expression modules. 862 863 Supplementary Tables 864 Supplementary Table S1 Population sampling, geographic coordinates, and experimental 865 design. 866 Supplementary S2 Basic information of the population’s vouchers, transcriptome sequencing, 867 and alignment rates of RNA libraries of four tissues of four different morphs of Stellera 868 chamaejasme. 869 Supplementary S3 Statistics of the Trinity assembly of transcriptomes for four sample tissues of 870 four different morphs of Stellera chamaejasme. 871 Supplementary Table S4 Statistics of the filtering of transcriptomes/pan - transcriptomes for 872 four sample tissues of four different morphs of Stellera chamaejasme. 873 Supplementary Table S5 OrthoFinder summary statistics for transcriptome assemblies across 874 floral morphs of Stellera chamaejasme. 875 Supplementary Table S6 Representative morph-associated DEGs and their functional roles in 876 pigment biosynthesis and regulatory processes. 877 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint 878 Figure 1. Eco -evolutionary framework linking floral polymorphism, regulatory 879 canalization, and hybrid predictability. 880 Conceptual framework illustrating how ecological divergence is associated with transcriptional 881 regulatory structure and expression outcomes in a mosaic morph (RWh). Floral morphs (PP, 882 RW, YW, YY) represent distinct phenotypic states hypothesized to correspond to canalized gene 883 regulatory networks, where co- expression modules exhibit reduced within- population variance. 884 A rare mosaic morph (RWh), occurring in zones of sympatry between PP and RW, represents a 885 perturbation of these regulatory states. Expression in RWh may be additive or parent -biased in 886 strongly canalized modules, whereas less constrained modules may exhibit non- additive or 887 transgressive expression. Thus, most modules remain stable, while a subset becomes 888 destabilized, generating transgressive regulatory states. Arrows indicate the progression from 889 ecological divergence to regulatory architecture and expression outcomes. 890 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint 891 Figure 2. Geographic distribution, phenotypic differentiation, and phylogenetic 892 relationships of floral morphs in Stellera chamaejasme. 893 (A) Geographic distribution of sampled populations across southwestern China. Sampling sites 894 are indicated by flower icons colored by dominant morph; background shading represents 895 elevation (m a.s.l.). Inset shows the study region within China . (B) Mean floral reflectance 896 spectra (300–700 nm) for four morphs (PP, RW, YW, YY), showing distinct spectral profiles 897 across the visible range . (C–G) Representative inflorescences of the five morphs: ( C) pure pink 898 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint (PP), (D) red–white (RW), (E) red–white mosaic (RWh), ( F) yellow–white (YW), and ( G) pure 899 yellow (YY). Symbols correspond to morph designations used throughout . ( H) Morph-level 900 phylogeny inferred from genome -wide single -copy orthologs (OrthoFinder), showing 901 relationships among morphs. The yellow lineage (YW, YY) diverges first, followed by 902 separation of RW from the pink lineage, with PP and RWh sharing the most recent common 903 ancestry. 904 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint 905 Figure 3. Gene expression divergence and functional differentiation among floral morphs. 906 (A) Principal component analysis (PCA) of normalized gene expression profiles showing 907 clustering by floral morph (PP, RW, RWh, YW, YY). Points represent individuals colored by 908 morph. PC1 and PC2 explain 21% and 11% of total variance, respectively, revealing strong 909 morph-associated transcriptional structure . (B) Venn diagram showing the number and overlap 910 of morph- associated differentially expressed genes (DEGs) among morphs. Numbers indicate 911 gene counts, with percentages relative to the total number of DEGs. Representative genes from 912 key metabolic and regulatory pathways are indicated. ( C) Volcano plot of differential gene 913 expression based on a generalized linear model (see Methods). Each point represents a gene 914 plotted by log₂ fold change and −log₁₀ adjusted p -value. Positive and negative log₂ fold change 915 values indicate relatively higher expression in one or the other morph depending on the model 916 contrast. Significantly up- regulated genes are shown in red, down- regulated genes in blue, and 917 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint non-significant genes in gray. Pairwise comparisons are shown in Supplementary Fig. S2. ( D) 918 KEGG pathway enrichment of DEGs. Dot size represents gene counts and color indicates 919 adjusted p- values. Enriched pathways include amino acid biosynthesis, carbon metabolism, 920 signal transduction, and secondary metabolism associated with pigmentation and stress 921 responses. 922 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint 923 Figure 4. Canalization and modular destabilization of gene regulatory networks across 924 floral morphs. 925 (A) Hierarchical clustering of genes based on expression similarity, with co- expression modules 926 identified using dynamic tree cutting. Colored bars indicate module assignments . (B) Canonical 927 analysis of principal coordinates (CAP) based on population- mean module eigengenes. Points 928 represent morph means (PP, RW, RWh, YW, YY), showing morph- associated structure in 929 regulatory profiles. ( C) Variance partitioning of module eigengene variation into components 930 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint associated with morph identity, geographic structure, their shared contribution, and residual 931 variance. (D) Within-morph expression variance (standard deviation, SD) across co -expression 932 modules (n = 53). Boxplots show reduced variance in parental morphs and elevated variance in 933 subsets of modules in RWh. Adjusted P -values (Benjamini–Hochberg) are indicated. (E) Hybrid 934 variance inflation across modules, calculated as the ratio of variance in RWh relative to the 935 parental mean. Modules are ordered by increasing inflation, showing that most remain stable 936 while a subset exhibit pronounced increases. ( F) Hybrid variance inflation across inheritance 937 modes (additive, dominant PP, dominant RW, transgressive). Points represent individual 938 modules, with higher variance inflation primarily associated with transgressive inheritance. 939 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint Table 1. Genome-wide inheritance patterns based on expression hybrid index. 940 Category Gene count Raw (n) Raw (%) Scaled (n) Scaled (%) Intermediate 4,446 8.93 3,463 6.89 PP-like 3,530 7.09 2,476 4.93 RW-like 4,172 8.38 1,085 2.16 Transgressive below PP 24,643 49.50 37,343 74.33 Transgressive above RW 12,995 26.10 5,871 11.69 Total transgressive 37,638 75.60 43,214 86.02 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint Table 2. Eco-evolutionary architecture of transcriptomic divergence, canalization, and RWh 941 expression predictability. 942 Analysis block Metric Statistic Value Interpretation Selection vs drift PERMANOVA (Morph) R² 0.45 Strong morph-structured divergence F 1.86 p 0.001* Highly significant dbRDA (spatial control) Overall model F 1.43 Morph effect after geography p 0.07 Marginal after spatial correction dbRDA (marginal) Morph F 1.43 p 0.08 Weak but persistent morph signal Variance partition Unique morph fraction (X1|X2) Adj. R² 0.16 Divergence consistent with selection Unique geography fraction (X2|X1) Adj. R² 0.12 Isolation-by-distance / drift Shared fraction Adj. R² 0.05 Spatially structured selection Hybrid predictability Additive modules % (n) 14 (26.42%) RWh expression approximates mid-parent values Dominant (PP) % (n) 15 (28.30%) Expression biased toward PP parent Dominant (RW) % (n) 4 (7.55%) Expression biased toward RW parent Transgressive modules % (n) 20 (37.74%) Novel regulatory states Hybrid canalization Inflation index (median) – 0.57 Most modules stable Inflation index (95th percentile) – 4.77 Subset highly unstable Canalization– inheritance coupling Inflation (Transgressive vs others) Wilcoxon p 0.25 No global coupling PERMANOVA and dbRDA based on population-mean module eigengenes (n = 14 populations). 943 Inflation index = SDRWh / mean parental SD(PP, RW, YW, YY). 944 Significance codes: *** P ≤ 0.001; · P ≤ 0.1. 945 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint Table 3. Top five transcriptional modules showing the strongest hybrid variance inflation and 946 distinct inheritance modes. 947 Module Inheritance mode SDPP SDRW SDRWh SDYW SDYY Parental mean SD Inflation index MEskyblue3 Transgressive 0.0273 0.0134 0.4730 0.0106 0.0087 0.0150 31.52 MEmagenta Transgressive 0.0351 0.0087 0.4420 0.0047 0.0147 0.0158 27.90 MEdarkolivegreen Transgressive 0.0411 0.0309 0.2310 0.0256 0.0259 0.0309 7.47 MEorange Dominant_PP 0.0365 0.0233 0.0680 0.0108 0.0206 0.0228 2.98 MEskyblue2 Additive 0.0545 0.0175 0.0820 0.0458 0.0341 0.0380 2.15 Hybrid variance inflation indices were calculated for the transgressive morph ( RWh) relative to 948 the mean within-population variance of the two closest parental morphs ( PP and RW). Modules 949 were classified as additive, dominant, or transgressive based on hybrid module eigengene means 950 relative to parental morphs. Standard Deviation (SD) values represent within-population standard 951 deviations of module eigengenes. 952 Inflation index = SD RWh / mean parental SD (PP, RW, YW, YY ); Parental morphs ( PP and RW) were 953 identified as the closest parental types to RWh based on module-level expression distance. 954 .CC-BY-NC 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted April 21, 2026. ; https://doi.org/10.64898/2026.04.16.719044doi: bioRxiv preprint