# Source paper archive - Source title: Heat Shock Transcription Factor B1 Coordinates Distinct bHLH Modules to Orchestrate Land Plant Adaptation - Author(s): Tien, H.-C., Lin, K.-H., Yang, H.-J., Yen, M.-R., Thalimaraw, L., Li, Y.-R., Cheng, C.-Y., Liu, M.-J., Wu, T.-Y. - Source landing URL: https://www.biorxiv.org/content/10.64898/2026.11.09.765476v1 - PDF URL: https://www.biorxiv.org/content/10.64898/2026.11.09.765476v1.full.pdf - PDF SHA-256: 694cc0f841f36da4e1245bf59482ca911f23615f1f426307015ed059a45db645 - Converter: pymupdf4llm 0.0.27 - Limits: Automated PDF text extraction; reading order, equations, tables and formatting may be incomplete or incorrect. Images are omitted, and no OCR was performed. Markup and links are neutralized for safe downloading. Consult the original PDF for the authoritative version. --- [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 1 **Title: Heat Shock Transcription Factor B1 Coordinates Distinct bHLH Modules to** 2 **Orchestrate Land Plant Adaptation** 3 Huan-Chi Tien [1,2,3#], Kuan-Hung Lin [1#], Lavakau Thalimaraw [1,4], Ming-Ren Yen [1], Hsiu-Ju 4 Yang [1], Ya-Ru Li [5], Chia-Yi Cheng [4,6], Ming-Jung Liu [5], and Ting-Ying Wu [1,2,3,4,7*] 5 1 Institute of Plant and Microbial Biology, Academia Sinica, Taipei 115201, Taiwan 6 2 Molecular and Biological Agricultural Sciences Program, Taiwan International Graduate 7 Program, National Chung Hsing University and Academia Sinica, Taipei 115201, Taiwan 8 3 Graduate Institute of Biotechnology, National Chung Hsing University, Taichung 402202, 9 Taiwan 10 4 Genome and Systems Biology Degree Program, Academia Sinica and National Taiwan 11 University, Taipei 10617, Taiwan 12 5 Biotechnology Center in Southern Taiwan, Academia Sinica, Tainan, 711010, Taiwan 13 6 Department of Life Science, National Taiwan University, Taipei 10617, Taiwan 14 7 Biotechnology Center, National Chung Hsing University, Taichung 402202, Taiwan 15 # These authors contribute equally to this study 16 *Correspondence: tingying@gate.sinica.edu.tw 17 18 Ting-Ying Wu 19 Institute of Plant and Microbial Biology 20 Academia Sinica 21 Taipei, Taiwan 22 tingying@gate.sinica.edu.tw 23 ORCID: 0000-0001-6488-5133 1 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 24 **Abstract** 25 Adaptation to terrestrial environments required the evolution of morphological structures, 26 such as rhizoids and air pores, and robust tolerance mechanisms to environmental stress. 27 While these innovations are well-described, the regulatory logic coordinating 28 development with stress signaling during this transition remains largely unknown. Here, 29 we identify Heat Shock Transcription Factor B1 (HSFB1) as a central evolutionary node 30 that bridges this gap in _Marchantia polymorpha_ . Using single-nucleus RNA sequencing 31 (snRNA-seq) and gene regulatory network analysis, we reveal that MpHSFB1 is required 32 to activate cell-type-specific programs for terrestrial adaptation, functioning via the 33 regulation of cytokinin and reactive oxygen species (ROS) homeostasis. We find that 34 MpHSFB1 functions through a modular logic: it co-opts a specific set of bHLH 35 transcription factors (including MpbHLH21/40) to drive developmental morphogenesis, 36 while recruiting a distinct partner (MpbHLH51) to confer resistance to prolonged heat 37 stress. This functional uncoupling allows a single ancestral regulator to manage disparate 38 physiological needs. By reconstructing HSFB1 regulons across diverse plant lineages, 39 we propose that MpHSFB1 acted as an ancestral pleiotropic hub that was effectively 40 rewired and subfunctionalized in later angiosperm lineages for land plant evolution. 41 **Keywords:** terrestrial adaptation, Marchantia polymorpha, cytokinin, snRNA-seq, 42 regulatory network evolution, prolonged heat stress 43 2 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 44 **Introduction** 45 Land plants underwent several key morphological modifications to adapt to terrestrial 46 environments. One such innovation was the development of rhizoids—root-hair-like 47 structures that arise from epidermal cells to facilitate water and nutrient absorption. 48 Another was the formation of air pore complexes on the plant surface, which enable gas 49 exchange while minimizing water loss [1] . In addition, land plants evolved mechanisms to 50 cope with oxidative stress caused by factors such as heat, water, nutrient scarcity, and 51 UV radiation. These mechanisms involve the perception of stressors and activation of 52 signal transduction pathways that include reactive oxygen species (ROS), specific 53 phytohormones, and the co-option of multiple transcription factors (TFs) [2] . The integration 54 of these processes requires the coordination of both internal and external signals [3] . 55 Consequently, understanding how plants perceive and process environmental cues 56 through multilayered molecular networks, and how these networks influence development 57 and survival strategies, remains a central focus of research on land plant adaptation. 58 One such group of phytohormones, cytokinin, plays a crucial role not only in plant 59 development but also in responses to abiotic and biotic stresses. Bryophytes already 60 possess a complete cytokinin signaling system, in which the signal is perceived and 61 transmitted through a variant of the Two-Component System (TCS), which consists of 62 CHASE domain–containing histidine kinase receptors (CHKs), histidine-containing 63 phosphotransfer proteins (HPTs), and response regulators (RRs) [4,5] . Although this 64 primitive signaling pathway originated in green algae, it acquired pleiotropic functions 65 during land colonization by integrating with other signaling and gene regulatory modules, 66 enabling the coordination of growth, development, and stress responses under terrestrial 67 conditions [6,7] . Consequently, mutants defective in cytokinin signaling exhibit diverse 68 morpho-physiological phenotypes across land plants, ranging from bryophytes to 69 angiosperms, including altered root and gametophyte development and increased 70 sensitivity to drought or salt stress [8] . However, the specific TFs that integrate these 71 developmental and environmental cues to facilitate adaptation remain largely unresolved. 72 While algal lineages possess ancestral Heat Shock Transcription Factor (HSF), the 73 distinct divergence into class A and class B HSFs first occurred in the ancestral land 3 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 74 plant _[9,10]_ . Unlike HSFA1, which serves as the conserved master regulator of high 75 temperature responses in plants, HSFB1 acquired more diversified functions [11-13] . For 76 example, Mp _hsfb1_ knockout mutants in _Marchantia_ exhibit only minor sensitivity to high 77 temperatures _[9]_ . In contrast, the development of air pores, air chambers, and rhizoids was 78 defective in _Marchantia_ Mp _hsfb1_ mutants, resembling the defects seen in mutants with 79 impaired cytokinin signaling [14-16] . Given the involvement of HSFB1 in traits related to 80 morphological adaptation to the terrestrial environment, we hypothesize that MpHSFB1 81 (Mp4g12230) evolved as a central regulator that balances external cues (e.g., 82 temperature) with internal hormonal signaling pathways (e.g., cytokinin). This 83 coordination likely orchestrates developmental responses and defense strategies 84 necessary to overcome the challenges of life on land. However, the underlying 85 mechanisms by which HSFB1 achieves such precise spatiotemporal regulation under 86 various conditions remain unclear. 87 Basic helix–loop–helix (bHLH) proteins are ancient TFs found across all eukaryotes and 88 represent one of the most expansive and evolutionarily versatile TF families in plants [17] . 89 Plant bHLHs form a monophyletic group distinct from those in rhodophytes (red algae) 90 and metazoans [18] . A major radiation of the bHLH family occurred after the divergence of 91 the moss and vascular plant lineages, giving rise to approximately 26–32 subfamilies that 92 largely emerged during the critical phase of land colonization [19] . This expansion likely 93 provided the regulatory versatility needed for plants to evolve novel organs and 94 developmental plasticity in fluctuating terrestrial environments. bHLHs commonly form 95 homo- or heterodimers and function as transcriptional activators and/or repressors of 96 genes involved in key developmental processes and their integration with environmental 97 cues such as high temperature and nutrient deficiency [17] . For example, the bHLHs ICE1/2, 98 SPCH, MUTE, and FAMA form heterodimers to regulate stomatal development, while 99 LRL1/2 and RSLs are involved in root hair formation in _Arabidopsis_ [20-22] . This 100 developmental role extends to bryophytes, where ICE1 (MpbHLH16, Mp4g04910) and 101 ZOU1 (MpbHLH38, Mp4g19650) promote pegged rhizoid differentiation, while RSL1 102 (MpbHLH14, Mp3g17930) drives general rhizoid formation [23,24] . In addition, illustrating 103 their role in environmental integration, the bHLH PIF4 regulates SPCH in response to 104 high ambient temperatures in _Arabidopsis_ [25] . Yet, how these bHLHs act synergistically or 4 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 105 antagonistically to regulate specific downstream targets, despite sharing similar DNA 106 binding sites, remains an open question. One plausible explanation is that they cooperate 107 with distinct TFs in specific regulatory contexts. Elucidating the combinatorial regulation 108 between bHLHs and other TFs will provide crucial insights into how the diversification of 109 transcriptional networks contributed to developmental complexity during land plant 110 evolution. 111 The ability to integrate these pathways is key to survival and colonization on land. In this 112 study, we demonstrate that HSFB1 functions as an ancestral regulatory hub that 113 maintains homeostasis between development and stress responses through modular 114 interactions with different bHLH TFs. Through integrated analyses of single-nucleus RNA 115 seq, bulk time-course RNA-seq, and phenotypic observations using various genetic 116 mutants, we found that MpHSFB1 and MpbHLHs regulate early developmental processes 117 in _Marchantia_ via cytokinin. In contrast, MpHSFB1 modulates a distinct set of MpbHLHs 118 to confer tolerance during prolonged heat stress. Furthermore, swapDAP-seq analysis 119 revealed that HSFB1 spatiotemporally regulates diverse gene functions across land 120 plants. Together, these findings highlight how HSFB-mediated diversification contributed 121 to the evolution of complex developmental and stress-response networks in plants. 122 **Results** 123 **MpHSFB1 regulates morphological adaptations associated with terrestrial** 124 **lifestyles in** _**Marchantia**_ 125 _Marchantia_ evolved various morphological and physiological adaptations to overcome 126 unfavorable environmental conditions. These structures included the development of air 127 pores, rhizoids, and mechanisms for hormone homeostasis and ROS regulation (Figure 128 1A). We observed developmental defects in the Mp _hsfb1_ knockout mutant during the 129 early developmental stage of _Marchantia_, starting from day 7 (Figure 1B, Supporting 130 Figure S1, Supporting Dataset S1). At day 14, the Mp _hsfb1_ mutant exhibited several 131 distinct morphological abnormalities compared to wild-type Tak-1 plants. Consistent with 132 previously reported phenotypes [9,26], these included a ~50% higher rhizoid-to-thallus 133 density, an absence of air pore and air chamber formation, and an abnormally upward 134 curved thallus. In addition to these known traits, we revealed several novel phenotypic 5 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 135 characteristics in the Mp _hsfb1_ mutant, including a ~3.5-fold smaller surface area, a ~1.3 136 fold lower solidity value, and ~1.8-fold higher ROS levels (Figures 1C and 1D, Supporting 137 Figure S1, Supporting Dataset S1). In contrast, the Mp _HSFB1_ ox transgenic lines showed 138 no differences in air pore and air chamber formation, thallus area, or rhizoid-to-thallus 139 ratio when compared to Tak-1 plants. However, they exhibited significantly higher solidity 140 (~1.2-fold) and lower ROS levels (~1.3-fold) (Figures 1C and 1D, Supporting Figure S1, 141 Supporting Dataset S1). Together with the evolution of HSFB1 in the ancestral land plant 142 lineage [9,10] (Supporting Figure S1C), we hypothesize that MpHSFB1 plays a crucial role 143 in regulating terrestrial adaptation. To test this, we designed a comprehensive study 144 integrating time-course RNA-seq, snRNA-seq, and DAP-seq analyses in Tak-1 (wild-type) 145 plants, as well as in Mp _hsfb1_ and Mp _HSFB1ox_ transgenic lines. These molecular data 146 were combined with time-series phenotypic observations to capture the regulatory roles 147 of MpHSFB1 across different developmental stages and cell types in _Marchantia_ (Figure 148 1E). 149 **MpHSFB1 influences cell type abundance and similarity in** _**Marchantia**_ 150 To investigate cellular-level changes in Mp _hsfb1_ mutants, we generated single-nucleus 151 RNA-seq (snRNA-seq) datasets from Tak-1 plants and Mp _hsfb1_ mutants using the 10x 152 Genomics platform (Supporting Figure S2A, Supporting Dataset S2). In total, we collected 153 over 20,000 nuclei per genotype. After quality control filtering, we retained 18,073 and 154 14,853 high-quality nuclei from Tak-1 and Mp _hsfb1_, respectively. The datasets showed 155 comparable depth, with a median of ~2,060 unique transcripts and ~1,150 detected genes 156 per nucleus for both genotypes (Figure 2A, Supporting Dataset S2). 157 We clustered these nuclei into distinct cell types based on a previously published 158 _Marchantia_ reference cell atlas [26] (Supporting Figure S2B, Supporting Dataset S2). Most 159 major cell types were successfully resolved and comparable to the protoplast-based 160 scRNA-seq atlas. We observed significant shifts in cell type abundance: rhizoid and 161 mucilage hair cells were significantly increased in the mutant, whereas gemma and 162 differentiating cells were markedly reduced (Figure 2B, Supporting Figure S3, Supporting 163 Dataset S2). This cellular composition shift mirrors the observed macroscopic phenotypes. 6 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 164 Furthermore, cell-type marker expression analysis revealed that air pore- and oil body 165 specific genes exhibited the greatest transcriptional divergence between Tak-1 and 166 Mp _hsfb1_, suggesting that cell-type identity is compromised in the mutant (Figure 2C, 167 Supporting Dataset S2). Consistent with this transcriptomic profile, histochemical staining 168 confirmed that Mp _hsfb1_ mutants possess a significantly reduced number of oil bodies 169 (Supporting Figure S3D & S3E). Pseudotime analysis indicated that the developmental 170 trajectories of rhizoid and differentiating cells were significantly altered in Mp _hsfb1_ 171 compared to Tak-1 (Figure 2D, Supporting Dataset S3). Finally, we examined the 172 differentially expressed transcription factor (TF) families between Tak-1 plants and 173 Mp _hsfb1_ mutants. We found that bHLH ( _n_ = 12, p < 0.01) and NAC ( _n_ = 3, p < 10 [-4] ) TFs 174 were among the most highly ranked families in both unique and common cell populations 175 (Figure 2E, Supporting Figures S3B & S4, Supporting Dataset S3). These results suggest 176 that MpHSFB1 regulates development in _Marchantia_ in a cell-type-specific manner, likely 177 through the co-regulation of bHLH and NAC transcription factors. 178 **The Mp** _**hsfb1**_ **regulons underwent functional rewiring through altered gene** 179 **targeting** 180 To further explore regulation at the cellular level, we constructed gene regulatory networks 181 (GRNs) to dissect the distinct and overlapping regulons in each genotype. Although ~75% 182 of the regulons ( _n_ = 32) were shared between the two genotypes, the networks from wild 183 type Tak-1 plants displayed a denser, more interconnected structure. Specifically, Tak-1 184 networks contained ~1.3-fold more nodes and exhibited ~1.5-fold higher average 185 connectivity compared to those in the Mp _hsfb1_ mutants (Figure 3A, Supporting Figures 186 5A-5C, Supporting Dataset S4). This suggests that the Mp _hsfb1_ mutant has a fragmented 187 regulatory system, leading to uncoordinated development and disrupted stress 188 homeostasis. 189 This structural difference was further supported by the observation that, while common 190 nodes accounted for 40% and 51% of the total nodes in Tak-1 and Mp _hsfb1_ regulons, 191 respectively, the Tak-1-specific regulons contained a significantly greater number of 192 nodes and edges (Figure 3B, Supporting Figure S5C, Supporting Dataset S4). 193 Additionally, the TF nodes in Tak-1 regulons displayed significantly higher betweenness, 7 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 194 Katz, and eigenvector centralities. This suggests that in the wild type, these TFs are highly 195 connected and play more integrative roles across the network compared to the mutant 196 (Figure 3C, Supporting Dataset S4). 197 We next examined the contribution of each regulon across cell types by ranking them 198 using the regulon specificity score (RSS). Among the 32 common regulons, 17 (53.1%) 199 retained their cell-type specificity, exhibiting the highest RSS rank in the same cell type in 200 both genotypes. However, 15 (46.9%) showed shifted specificity, ranking highest in 201 different cell types between the two genotypes (Figure 3D, Supporting Figure S5B, 202 Supporting Dataset S4). 203 Gene Ontology (GO) analysis further revealed that the biological functions of these 204 regulons diverged significantly. The regulons in Mp _hsfb1_ mutants were primarily 205 associated with housekeeping and structural functions such as photosynthesis, 206 translation, and cell wall organization. In contrast, regulons in Tak-1 plants were enriched 207 for functions linked to environmental adaptation, including cytokinin signaling and 208 oxidative stress responses (Figure 3E, Supporting Figure S3C, Supporting Dataset S5). 209 These results suggest that while the core regulatory architecture is partially conserved, a 210 substantial fraction of regulons has been functionally reassigned in the absence of 211 MpHSFB1, indicating HSFB1-dependent rewiring of cell-type-specific regulatory 212 programs. 213 **MpHSFB1 and MpbHLHs regulate cytokinin homeostasis during the early** 214 _**Marchantia**_ **development** 215 Developmental bulk RNA-seq at days 3, 7, and 14 identified 4,011 differentially expressed 216 genes (DEGs) across eight clusters (C1-C8) (Figure 4A, Supporting Dataset S6). In 217 Mp _HSFB1ox_ transgenic lines, cluster C7 genes peaked at day 7 and were significantly 218 enriched for cytokinin response functions (q-value < 10⁻⁸) (Figure 4B, Supporting Figures 219 S6A & S6B). Consistent with this, cytokinin biosynthesis genes were differentially 220 expressed between Tak-1 and Mp _hsfb1_ mutants across the developmental time course 221 and within specific tissues, including air pore cells, gemma cells, and ventral scales from 222 the snRNA-seq analysis (Supporting Figures S2-S4 & S6C, Supporting Dataset S2 & S6). 223 To determine whether this transcriptional signature was associated with altered cytokinin 8 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 224 homeostasis, we measured endogenous active cytokinin ( _cis_ -zeatin) levels. _Cis_ -zeatin 225 concentrations were significantly lower (~2-fold, _p_ = 0.001) in Mp _hsfb1_ mutants compared 226 to Tak-1, but approximately 2- to 7-fold higher in Mp _HSFB1ox_ transgenic lines across all 227 time points (Figure 4C, _p_ = 0.007). Furthermore, exogenous cytokinin application rescued 228 the Mp _hsfb1_ developmental phenotypes, restoring both thallus area (~2-fold increase) 229 and thallus angle (~1.3-fold increase) to near-wild-type levels (Figure 4D & 4E, Supporting 230 Dataset S1). Together, these results indicate that MpHSFB1 regulates early _Marchantia_ 231 development, at least in part, through cytokinin homeostasis. 232 We next sought to identify downstream transcriptional regulators that may connect 233 MpHSFB1 to these cytokinin-associated developmental responses. The bHLH TF family 234 was strongly enriched in the cytokinin-related clusters C6 and C7 (Figure 4F), while 235 MpHSFB1-co-expressed MpbHLHs were also enriched in rhizoid cells based on snRNA 236 seq (Fisher's exact test, _p_ < 0.05; Supporting Dataset S2). We therefore focused on the 237 regulatory relationship between MpHSFB1 and these MpbHLHs by integrating multiple 238 genomic and transcriptomic datasets. DAP-seq identified several established 239 developmental MpbHLHs, including MpbHLH4, 6, 16, and 37, as direct MpHSFB1-bound 240 targets (Figure 4G, solid arrows) [23,24,27,28] . In parallel, DEGs identified in public Mp _bhlh_ 241 mutant datasets significantly overlapped with those altered in our Mp _hsfb1_ mutant dataset, 242 including MpBHLH7, 21, and 40 (Supporting Dataset S6 & S7). MpHSFB1 expression, 243 however, remained unchanged in these Mp _bhlh_ mutant backgrounds, supporting a 244 regulatory direction from MpHSFB1 toward these bHLHs rather than the reverse. In 245 addition, snRNA-seq identified MpbHLH28 as a key regulon component in mucilage hairs, 246 which are structures involved in preventing water loss [29] (Figures 3D & 4G). Based on 247 these complementary analyses, we selected MpbHLH21, 28, and 40, which showed the 248 strongest combined evidence across these analyses, for downstream functional 249 characterization. 250 Luciferase reporter assays showed that MpHSFB1 was sufficient to repress the 251 MpBHLH21 and MpBHLH28 promoters, but not the MpBHLH40 promoter (Supporting 252 Figures S7 & S8, Supporting Dataset S6 & S7). These results indicate that MpHSFB1 253 regulates individual MpbHLHs through distinct mechanisms, with MpBHLH21 and 254 MpBHLH28 showing direct evidence of transcriptional regulation by MpHSFB1, whereas 9 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 255 the relationship with MpbHLH40 may involve additional regulatory components. 256 Furthermore, we overexpressed MpHSFB1 in Mp _bhlh_ knockout backgrounds to 257 genetically test whether these bHLH factors function downstream of MpHSFB1. We 258 compared developmental phenotypes of Tak-1, Mp _hsfb1-25_, Mp _HSFB1ox-1_, Mp _bhlh40-_ 259 _5_, Mp _bhlh28-18_, Mp _bhlh21-4_, and the corresponding double lines 260 (Mp _bhlh40_ /Mp _HSFB1ox_, Mp _bhlh28_ /Mp _HSFB1ox_, and Mp _bhlh21_ /Mp _HSFB1ox_ ). In all 261 three genetic backgrounds, MpHSFB1 overexpression failed to override the respective 262 Mp _bhlh_ mutant phenotypes: the double lines retained the phenotypes of the 263 corresponding Mp _bhlh_ mutants in thallus size, solidity, and rhizoid-to-thallus ratio 264 (Supporting Figure S9). The epistasis of the Mp _bhlh_ mutations over MpHSFB1 265 overexpression provides genetic evidence that MpBHLH21, MpBHLH28, and MpBHLH40 266 function downstream of MpHSFB1 in developmental regulation. 267 We then asked whether these downstream MpbHLHs contribute to the developmental 268 and cytokinin-related phenotypes associated with MpHSFB1. Mpb _hlh21_ mutants 269 displayed a significantly lower rhizoid-to-thallus ratio (~0.67-fold), whereas Mp _bhlh40_ 270 mutants exhibited a higher ratio compared to Tak-1 (Figure 5A & 5B, ~2.7-fold), 271 resembling the rhizoid phenotype of Mp _hsfb1_ . Beyond rhizoid development, Mp _bhlh28_ 272 mutants displayed deeply embedded gemmae within the gemma cups, whereas 273 Mp _bhlh40_ mutants failed to maintain dormancy (Supporting Figure S7B). Crucially, 274 cytokinin levels were inversely correlated with rhizoid density ( _ρ_ = -0.67): Mp _bhlh21,_ 275 which produced fewer rhizoids, accumulated significantly higher cytokinin levels, whereas 276 Mp _bhlh40,_ which produced more rhizoids, showed significantly lower cytokinin levels 277 (Figure 5C & 5D). Exogenous cytokinin partially rescued the Mp _bhlh40_ mutant phenotype 278 by restoring thallus area but had no detectable effect on Mp _bhlh21_ mutants (Figure 5E & 279 5F), further linking MpBHLH40 to the low-cytokinin, high-rhizoid phenotype observed in 280 the Mp _hsfb1_ mutant. 281 Consistent with these phenotypic differences, genotype-specific DEG analysis revealed 282 distinct regulatory functions for these two MpbHLHs. Mp _bhlh40_ DEGs were enriched in 283 cytokinin responsiveness, pattern formation, and root development, whereas Mp _bhlh21_ 284 DEGs were associated with sucrose metabolism, abiotic stress, and ROS pathways 285 (Figure 5G & 5H, Supporting Dataset S7). Integration with the Mp _hsfb1_ and Mp _HSFB1ox_ 10 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 286 time-course profiles further identified shared target genes within clusters C1, C4, and C6 287 that were misregulated in both lines and enriched in ROS regulation, lipid metabolism, 288 and hormone signaling (cytokinin and GA) (Supporting Dataset S7). Collectively, these 289 regulatory, genetic, and physiological data support a model in which Mp _BHLH21_ and 290 Mp _BHLH40_ function downstream of MpHSFB1 to coordinate distinct cytokinin-associated 291 developmental and physiological programs. 292 **MpHSFB1 regulates prolonged HS tolerance via modulating ROS response** 293 Regulating ROS homeostasis is essential for thermotolerance under prolonged HS [30] . 294 Under prolonged HS, Mp _HSFB1ox_ lines exhibited accelerated pigmentation and 295 significantly reduced thallus area compared to Tak-1 and Mp _hsfb1_ plants (Figure 6A, 296 Supporting Figure S10). While ROS levels decreased significantly in Tak-1 and Mp _hsfb1_ 297 lines post-HS ( _p_ = 1.8 × 10⁻¹³ and 0.018, respectively), they remained stable in 298 Mp _HSFB1ox_ lines, suggesting that overexpression prevents proper ROS suppression 299 and triggers oxidative stress under prolonged HS (Figure 6C, Supporting Figures S10C 300 & S10D, Supporting Dataset S1). RNA-seq identified hundreds of genotype-specific 301 DEGs, including 35 with opposite expression patterns between genotypes (Figure 6D, 302 Supporting Figure S10B, Supporting Dataset S8). These targets were enriched in ROS 303 responses, oxidative stress, and pigment biosynthesis (Figure 6E), supporting the 304 hypothesis that MpHSFB1 plays a pivotal role in maintaining ROS homeostasis during 305 temperature adaptation. 306 Expression profiling revealed that Mp _BHLH51_ and Mp _NAC6_ (MpJINGASA [31] ) exhibited 307 the most pronounced differences in expression between Mp _HSFB1ox_ lines and Mp _hsfb1_ 308 mutants (Figure 6F, Supporting Figure S10E). Functional analyses of their respective 309 mutants revealed divergent regulatory roles, where Mp _bhlh51_ mutants developed 310 significantly larger thalli (Figures 6G & 6H) and accumulated less ROS under stress 311 (Figure 6I, Supporting Figure S10F), consistent with an enrichment of their DEGs in 312 generalized stress adaptations such as ROS detoxification, cell wall reorganization, and 313 phenylpropanoid metabolism (Figures 6J & 6K, Supporting Dataset S8). Conversely, 314 Mp _nac6_ mutants displayed severe hypersensitivity to prolonged HS, characterized by 315 restricted growth and elevated ROS accumulation, which aligns with the enriched function 11 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 316 that was specialized in redox homeostasis and osmotic stress (Figure 6G-6K). Together, 317 these contrasting downstream cascades demonstrate that MpHSFB1 acts as a central 318 coordinator of thermal adaptation by fine-tuning distinct, specialized regulatory modules. 319 To test the regulatory specificity of downstream pathways, we performed cross 320 phenotypic evaluations of the development- and stress-related mutants. When 321 challenged with prolonged HS, the development-associated Mp _bhlh21_ and Mp _bhlh40_ 322 mutants maintained a thallus area and ROS ratio indistinguishable from wild-type Tak-1 323 plants, indicating they did not participate in prolonged HS (Supporting Figure S11A & 324 S11B). Reciprocally, under control conditions, the stress-responsive Mp _bhlh51_ and 325 Mp _nac6_ mutants, as well as Mp _bhlh51_ /Mp _HSFB1ox_ lines, exhibited no significant 326 abnormalities in their rhizoid-to-thallus ratio or endogenous cytokinin concentrations 327 (Supporting Figure 9, Supporting Figure S11C - S11E). Together, these results support 328 distinct functional roles for these downstream regulators, with MpBHLH21 and 329 MpBHLH40 primarily associated with developmental and cytokinin-related processes, 330 and MpBHLH51 and MpNAC6 with prolonged HS responses. 331 **MpHSFB1 and its orthologs coordinate tissue- and cell type-specific transcriptional** 332 **programs in plants** 333 Co-expression networks across species revealed that distinct HSFBs serve as regulatory 334 hubs within different cell types (Supporting Dataset S9). For example, HSFB2a was top 335 ranked in _Arabidopsis_ meristem procambium cells and the rice epidermis, while HSFB1 336 acted as a hub in tomato and maize root cap cells (Supporting Dataset S9). Based on GO 337 term similarity, the functions of the genes associated with each HSFB module were 338 grouped into seven clusters (Supporting Figure S12A, Supporting Dataset S9). Focusing 339 on these functions, MpHSFB1 showed a broad, even distribution across GO groups, 340 whereas angiosperm HSFB modules displayed narrower, more specialized functional 341 distributions. Specifically, genes associated with AtHSFB2a, SlHSFB1, and SlHSFB3a 342 exhibited a significantly higher proportion of stress-related functions and a depletion of 343 development-associated functions ( _p_ < 0.05 compared to MpHSFB1) (Supporting Figure 344 S12B). These data suggest that the ancestral MpHSFB1 regulates a broad, evenly 12 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 345 distributed array of functional target pathways, whereas angiosperm HSFBs have shifted 346 toward more specialized functional target profiles. 347 We next investigated the conservation of HSFB regulatory logic at the cellular level. Since 348 the expression patterns of a TF’s conserved target genes allow us to chart TF activity 349 across species [32], we performed swapDAP-seq to assess the binding capacity of 350 MpHSFB1 on the _Arabidopsis_ genome (Figure 7A, Supporting Figure S13A, Supporting 351 Dataset S10). By profiling MpHSFB1 binding to both Marchantia and _Arabidopsis_ 352 genomic DNA, we identified 535 target genes (392 orthogroups) and 2,302 target genes 353 (2,082 orthogroups), respectively (Figure 7B). Only 63 overlapping orthogroups (2.6%) 354 were shared between the two species, which were strictly enriched for core heat stress 355 response (FDR = 1.85 × 10 [-6] ) and ROS response (FDR = 0.011). Conversely, the closest 356 _Arabidopsis_ homolog, AtHSFB4, bound far fewer genomic sites, yielding no overlapping 357 orthogroups between genomes (Figure 7C) [33] . Furthermore, its _Arabidopsis_ -specific 358 targets were tightly restricted to specific localized functions, such as malate transport 359 (FDR = 5.93 × 10 [-3] ) and cell division (FDR = 0.013) (Figure 7D). Thus, genomic targeting 360 has diverged within the HSFB clade, where MpHSFB1 retains broad developmental and 361 stress-responsive binding potential, whereas tested angiosperm HSFBs display 362 significantly narrower profiles. 363 Finally, we validated these regulatory shifts using publicly available snRNA-seq and 364 scRNA-seq datasets alongside our own _Marchantia_ snRNA-seq data (Supporting Figure 365 S14, Supporting Dataset S10). While Arabidopsis HSFB regulon activity is broadly 366 enriched in undifferentiated meristematic cells across seedling, leaf, and root tissues, 367 MpHSFB1 activity is strongly associated with specialized, differentiated cell types like air 368 pores and oil body cells (Figure 2). Notably, MpHSFB1 shared regulatory features with 369 both stress-responsive cell types (resembling AtHSFB2a and AtHSFAs) and developing 370 cell types (resembling AtHSFB3 and AtHSFB4) (Supporting Figure S14). This cell-type 371 dependent versatility is associated with distinct HSFB–bHLH co-regulatory arrangements 372 (Figure 7E). In Arabidopsis roots, AtHSFB3 and AtHSFB4 correlated positively with cell 373 division-regulating TFs (AtbHLH074, AtbHLH77, and AtbHLH104). Similarly, MpHSFB1 374 exhibited dual regulatory associations: it correlated strongly with stress-related bHLHs 375 and HSFBs in seedling cells, but with development-related TFs in roots and leaves 13 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 376 (Figure 7E). These findings suggest a model where cell-type-specific HSFB–bHLH 377 interactions potentially underlie distinct adaptive traits in _Marchantia_ and _Arabidopsis_ 378 (Figure 7F). 379 **Discussion** 380 Overcoming terrestrial challenges was a pivotal step in land plant evolution, coinciding 381 with the diversification and expansion of the HSFB and bHLH transcription factor families. 382 These events enabled the regulatory innovations required for life on land. In this study, 383 we demonstrate that MpHSFB1 functions as an ancestral regulatory hub that coordinates 384 distinct bHLH and NAC transcription factor modules to orchestrate developmental and 385 stress-responsive programs across different cell types (Figures 2-7). By integrating 386 cytokinin signaling, MpHSFB1 modulates bHLH21 and bHLH40 to fine-tune the 387 spatiotemporal dynamics of rhizoid development, while independently coordinating 388 bHLH51 and NAC6 to mediate adaptive tolerance under prolonged heat stress. Such 389 modularization illustrates how conserved regulatory networks were progressively co 390 opted and rewired to generate new adaptive traits during land plant evolution. Collectively, 391 our findings suggest that ancestral HSFB transcriptional regulons were repurposed 392 through regulatory rewiring to facilitate the specialized cellular functions during land plant 393 evolution. 394 The transition of plants to land required the emergence of new cell types to form tissues 395 and organs, accompanied by the evolution of sophisticated tissue-to-tissue 396 communication systems. The evolution of GRNs involving cytokinin signaling likely 397 facilitated the coordinated regulation of stress acclimation and developmental processes 398 in land plants (Figures 3, 5, and 7) [34], including the plastic development of root or rhizoid 399 systems and the accumulation of protective pigments such as anthocyanins or 400 auronidins [35] . Indeed, endogenous cytokinin concentrations were misregulated in 401 Mp _hsfb1_ and Mp _bhlh_ mutants, and exogenous cytokinin application partially rescued their 402 developmental phenotypes (Figure 5), supporting the hormone's role in integrating the 403 molecular cues underlying these adaptive traits. Furthermore, HS is known to reduce 404 cytokinin content, whereas exogenous cytokinin generally enhances HS resistance [36-38] . 405 Cytokinin and HS responses cross-talk by regulating HSPs and photosynthesis-related 14 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 406 proteins, thereby influencing stomatal movement [38] . Given the misregulation of cytokinin 407 signaling and cytokinin-responsive genes in Mp _hsfb1_ mutants, it is likely that the 408 integration of HSFB into stress- and hormone-responsive pathways functioned as a 409 critical adaptive signaling mechanism during plant evolution (Figure 7). 410 We propose that MpHSFB1 represents an evolutionarily ancient HSFB that originally 411 possessed multimodal functions, which later diverged through subfunctionalization within 412 the angiosperm lineage. In _Marchantia_, MpHSFB1 regulates both developmental and 413 stress-related processes (Figures 5 and 6). By contrast, the _Arabidopsis_ genome encodes 414 five HSFB genes with specialized roles: AtHSFB1, AtHSFB2a, and AtHSFB2b are 415 primarily involved in stress regulation, whereas AtHSFB3 and AtHSFB4 function mainly 416 in developmental control. For example, AtHSFB4 (SCHIZORIZA) regulates root cell-fate 417 segregation, while the closest homolog of AtHSFB3 in black cottonwood controls xylem 418 cell-wall biosynthesis [13,39] . Interestingly, despite sharing nearly identical binding motifs, 419 AtHSFB3 and AtHSFB4 targeted fewer genes in DAP-seq datasets than AtHSFAs or 420 AtHSFB2a. A similar trend was observed when comparing MpHSFA1 and MpHSFB1. In 421 our cross-species comparison, we found that MpHSFB1 is capable of binding to 422 _Arabidopsis_ genomic DNA, whereas AtHSFB4 cannot bind to the _Marchantia_ genome. In 423 contrast, HSFA1 proteins exhibit conserved binding capacities across species boundaries 424 (Figure 7, Supporting Figure S13). Consistent with this, the HSFB lineage underwent 425 strong positive selection in bryophytes relative to the ancestral HSFs found in algae, likely 426 driven by the drastic environmental shifts associated with land plant adaptation [10] . These 427 results indicate that sequence divergence within the HSFB DNA-binding domain (DBD) 428 was a key determinant of the functional diversification and specialization observed during 429 plant evolution. 430 While direct protein–protein interactions between MpHSFB1 and the MpbHLHs were 431 undetected (Supporting Figure S7), DAP-seq revealed that the identified motifs were 432 generally longer in Marchantia (Figure 5) _[40,41]_ . Because many TF families form higher 433 order complexes to recognize longer motifs and increase DNA-binding specificity, these 434 longer motifs suggest the HSFB family may form similar complexes to enhance targeting 435 specificity and enable modular control [42-44] . Coupled with evidence indicating that HSFBs 15 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 436 cooperate with distinct bHLH TFs in a cell-type-specific manner (Figure 7), these findings 437 illustrate how TFs can be co-opted to generate new, cell-type-specific regulatory functions. 438 Like other eukaryotic TFs, bHLH factors form distinct regulatory complexes to govern 439 diverse developmental and stress responses. For example, the _Arabidopsis_ homologs of 440 MpBHLH40 (AtbHLH037 and AtbHLH088) control reproductive development and 441 photomorphogenesis [45], while MpBHLH21 and MpBHLH28 homologs regulate cell 442 elongation and root hair growth, indicating this functional diversity had already emerged 443 in bryophytes. This regulatory versatility may be driven by several mechanisms: 444 association with ACT-like domains, post-translational control, or the cell-type-specific 445 modulation of hormone pathways like cytokinin signaling [46] . Together, these features 446 underscore how bHLH TFs drove the diversification of plant developmental and stress 447 responsive programs. In conclusion, these results suggest how ancestral HSFB-mediated 448 regulatory modules evolved to facilitate cell-type specificity and functional diversity in land 449 plants. 450 **Material and Methods** 451 **Phylogenetic tree construction** 452 Full-length protein sequences were retrieved from the Tapscan database [47], JGI 453 PhycoCosm [48], and recent genome releases [49] . Our dataset encompasses key 454 evolutionary lineages, ranging from streptophyte algae ( _Klebsormidium nitens [v1.0]_, 455 _Mesotaenium endlicherianum_ _[SAG12.97]_, _Zygnema circumcarinatum [Zygcir6981b_2]_ ) 456 to bryophytes ( _Anthoceros agrestis [v1.0]_, _Marchantia polymorpha [v6.1r2]_ ), 457 tracheophytes ( _Selaginella moellendorffii [v1.0]_, _Azolla filiculoides [v1.2]_, _Salvinia_ 458 _cucullata [v2]_, _Picea abies [v1]_, _Ginkgo biloba [v1]_, _Oryza sativa [MSU v7.0]_, _Zea mays_ 459 _[B73_v4]_, _Arabidopsis thaliana [Araport11]_, _Solanum lycopersicum [ITAG4.1]_ ), with 460 _Chlorophyta_ and _Rhodophyta_ species serving as outgroups. Conserved DBDs were 461 identified using the NCBI Conserved Domain Database. We performed multiple sequence 462 alignment via MAFFT (v7.520) using the L-INS-i algorithm (--maxiterate 1000 --localpair) 463 and reconstructed the maximum-likelihood phylogenetic tree using IQ-TREE (v2.2.2.6). 464 The optimal amino acid substitution model was selected via ModelFinder Plus based on 16 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 465 the Bayesian Information Criterion, with branch support evaluated using 1,000 Ultrafast 466 Bootstrap replicates. 467 **Growth condition and morphological characterization** 468 _Marchantia polymorpha_ wild type (accession Tak1), MpHSFB1 overexpression lines 469 (Mp _HSFB1ox_ ), MpHSFB1 knockout lines (Mp _hsfb1_ ), and TF candidate knockout mutants 470 were grown on 1/2 × Gamborg B5 medium (pH 5.5) with 1% agar. Plants were grown 471 under continuous light at either 22 °C or 30 °C, or with and without the treatment of 10 472 μM 6-Benzylaminopurine (BA, Sigma). Morphological characterization was performed in 473 two phases: developmental and treatment. For the developmental phase, plants were 474 cultured at 22 °C, and images were captured on days 3, 7, and 14. For the treatment 475 phase, plants were cultured at either 22 °C or 30 °C, with or without BA treatment, and 476 images were captured on day 14. All images were obtained using a Zeiss Lumar V12 477 microscope, while other images were captured using a digital camera and further 478 analyzed using Fiji ImageJ software. 479 **Phenotypic analysis** 480 _Thallus area and shape analysis:_ The scale of each image was calibrated by drawing a 481 straight line along the scale bar (Set Scale). Measurement parameters, including area, 482 shape descriptors, and area fraction, were selected. Each image was then split into red, 483 green, and blue channels (Split Channels), and subsequent analyses were performed on 484 the blue channel. To reduce image noise, a Gaussian blur filter was applied (Gaussian 485 Blur). Images were automatically thresholded to distinguish the thallus from the 486 background. Non-thallus regions were cropped to retain only the area of interest, and the 487 image was inverted to prepare for particle analysis. Thallus area and shape 488 characteristics were quantified using “Analyze Particles,” and the corresponding values 489 were recorded. The convex hull of each thallus area was generated (Convex Hull), and 490 measurements were repeated to obtain convex hull area values. 491 _Rhizoid area analysis:_ For rhizoid quantification, the scale was calibrated, and an area of 492 interest was selected. Each image was split into red, green, and blue channels, and either 493 the red or blue channel was used, depending on which provided better contrast between 494 rhizoids and the background. Thresholding and selection of rhizoid areas were performed 17 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 495 manually. Non-rhizoid regions were cropped to retain only the rhizoid area, and the image 496 was inverted prior to measurement. 497 _Angle measurement:_ Side-view images were used to measure plant shoot angles. Each 498 image was cropped to include only the green portion of the plant. The Angle tool in ImageJ 499 was used for measurement. For the left-side angle, three points were selected 500 sequentially: (1) the top-left corner of the plant, (2) the midpoint of the plant, and (3) a 501 point forming a horizontal line with the midpoint on the left side, after which the angle was 502 measured. For the right-side angle, the midpoint was maintained, while the first point was 503 moved to the top-right corner and the third point to the right side, again forming a 504 horizontal line with the midpoint. Measurements were repeated as described above. 505 **Overexpression and knockout lines in** _**Marchantia**_ 506 The MpHSFB1 overexpression construct was generated using Mol-Clo–based Golden 507 Gate assembly. The coding sequence (CDS) of MpHSFB1 (synthesized by 508 OmicBioscience) was fused to an N-terminal 3×Flag tag and a C-terminal mCherry and 509 cloned under the control of the CaMV35S promoter with a NOS terminator to generate a 510 Level 1 (Lv1) plasmid. This Lv1 plasmid was subsequently assembled with a second Lv1 511 plasmid containing a hygromycin resistance cassette driven by the MpEF1α promoter, 512 yielding the Level 2 (Lv2) plasmid used for Agrobacterium-mediated transformation [9] . To 513 complement Mp _bhlh_ mutants with MpHSFB1, the CmmT control construct 514 (CaMV35S:mCherry:YFP:miniTurbo:NOS) and the CBmT construct expressing 515 MpHSFB1 (CaMV35S:mCherry:MpHSFB1:miniTurbo:NOS) were generated via Golden 516 Gate assembly. S uccessful transformants were screened on 1/2 × B5 medium plates 517 containing 100 µg/mL cefotaxime or 10 mg/L hygromycin. Two independent 518 overexpression lines at the T2 generation were selected based on transgene expression 519 levels determined by qPCR. 520 The Mp _hsfb1_ mutant was generated in our previous publication [9] . Knockout mutants of 521 MpbHLH21, MpbHLH28, MpbHLH40, and MpbHLH51 were produced using the 522 CRISPR/Cas9 (clustered regularly interspaced short palindromic repeats/CRISPR 523 associated protein 9) genome-editing system. Two guide RNAs (gRNAs) targeting 524 sequences flanking a short DNA fragment within each locus were designed and cloned 18 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 525 into either the pMpGE010 (Hygromycin resistance) or pMpGE011 (chlorsulfuron 526 resistance) expression vector. Successful transformants were screened on 1/2 × B5 527 medium plates containing 100 µg/mL cefotaxime and 0.5 µM chlorosulfuron or 10 mg/L 528 hygromycin. The mutation sites of two independent knockout lines at the T2 generation 529 were confirmed via PCR and sequencing. Gene locus identifiers and primers used for 530 sequencing and cloning are listed in Supporting dataset S11. 531 **Detection of reactive oxygen species (ROS)** 532 The histochemistry staining was followed by a previous study [16] . _Marchantia_ thalli were 533 harvested from 14-day-old Tak-1, Mp _HSFB1ox_ lines, Mp _hsfb1_ lines, and TF candidate 534 knockout mutants and stained with 3,3ʹ-Diaminobenzidine (DAB). The 1 mg ml−1 DAB 535 staining solution was directly applied on the plate, followed by 5 min vacuum and shake 536 for 1 h, at 50rpm in the dark incubation. The DAB staining solution was replaced with 537 distilled water, subsequently with 75% ethanol to remove the chlorophyll. Plant images 538 were taken and analyzed by Fiji ImageJ software. 539 **Endogenous cytokinin detection** 540 Thalli from 3-, 7-, and 14-day-old Tak-1, Mp _HSFB1ox_, and Mp _hsfb1_ lines were ground 541 using a pestle and resuspended in Bieleski buffer (60% methanol, 25% CHCl 3, 542 10%HCOOH, and 5% H2O). The mixture was sonicated on ice for 30 min, followed by 543 centrifugation at 13,000 rpm for 10 minutes at 4°C. The supernatant was transferred into 544 a new 1.5 mL tube and dried using a vacuum concentrator at 4 °C. UPLC-MS/MS analysis 545 was performed by ABRC Metabolomics Core. Deuterium-labeled trans-zeatin (tZ-d5) was 546 used as an internal control. Absolute quantification of cytokinins was calculated using the 547 following formula: [cZ (area) /d5-tZ (area) x_0.5_ng (d5-tZ) ]÷FW (mg) =____ng/mg. 548 **Oil Body Staining** 549 For oil body observation, gemmae were collected from Tak-1 plants, Mp _HSFB1ox_ lines, 550 and Mp _hsfb1_ mutants. The fluorophore BODIPY 493/503 (4,4-difluoro-1,3,5,7,8 551 pentamethyl-4-bora-3a,4a-diaza-s-indacene; Thermo Fisher Scientific) was dissolved in 552 water to prepare a 0.5 mg/mL stock solution, which was stored at -20 °C. Prior to imaging, 553 the tissues were incubated in a 0.05 mg/mL BODIPY working solution for 10 min, followed 19 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 554 by two washes with ddH 2 O. Fluorescence was visualized using a Zeiss LSM 880 confocal 555 microscope with a 488 nm excitation argon laser. Finally, the number of oil bodies per 556 gemma was quantified manually. 557 **Protein-protein interaction (PPI) assay using RUBY-based reporter system** 558 The MpHSFB CDS was fused to the N-terminal GAL4 DNA-binding domain (GAL4) under 559 the control of a double CaMV35S promoter. Candidate TF CDSs were synthesized (by 560 Twist Bioscience), cloned in-frame with the C-terminal GAL4 activation domain (p65), and 561 driven by the double CaMV35S promoter. The three constructs were introduced into 562 Agrobacterium tumefaciens by electroporation. Agroinfiltration was performed as 563 previously described with minor modifications [50] . Briefly, transformed Agrobacterium 564 strains were cultured in 5 mL LB broth supplemented with 250 µg/mL carbenicillin, 10 mM 565 MES (pH 5.7), and 20 μM acetosyringone at 28 °C for 2 days. Cells were collected by 566 centrifugation and resuspended in infiltration buffer (10 mM MgCl₂, 10 mM MES, pH 5.7, 567 and 200 μM acetosyringone) to an OD₆₀₀ of 1.0, then incubated at room temperature for 568 at least 1 h. Equal volumes of the three Agrobacterium suspensions were mixed (1:1:1) 569 and infiltrated into fully expanded leaves of 4-week-old _Nicotiana benthamiana_ plants 570 using a needleless syringe. Betalain pigment accumulation was monitored 2 days post 571 infiltration. 572 **Generation and analysis of the DAP-seq dataset** 573 _Plasmid construction:_ T he coding sequences (CDS) of MpHSFB1 were synthesized by 574 Twist Bioscience, and the CDS of Arabidopsis HSFB4 (AtHSFB4) was amplified from 575 _Arabidopsis_ cDNA and subcloned to Lv0 acceptor pAGM9121. The CDSs were cloned 576 downstream of the SP6/T7 promoter with an N-terminal HaloTag and a T7 terminator to 577 generate level 1 (Lv1) expression vectors. Each CDS was first subcloned into pAGM9121 578 as a level 0 module using BbsI, and subsequently assembled into the pICH47732 Lv1 579 vector via Golden Gate cloning, yielding an SP6/T7 promoter–driven construct encoding 580 a HaloTag-fused CDS followed by a T7 terminator. Recombinant proteins were expressed 581 in vitro using the TNT SP6/T7 High-Yield Wheat Germ Protein Expression System 582 (Promega) according to the manufacturer’s instructions. Protein expression and purity 583 were verified by SDS–PAGE. 20 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 584 _DNA Fragmentation and Preparation:_ Approximately 100 mg of _Marchantia_ gemmaling 585 (14 days old) were harvested and ground in liquid nitrogen. Genomic DNA was extracted 586 using the CTAB-phenol/chloroform method and sheared mechanically to an average size 587 of 300–500 bp. DNA fragments were end-repaired, A-tailed, and ligated to indexed 588 sequencing adapters for multiplexing. Library preparation, including fragmentation and 589 adapter ligation, was performed by the Genomic Technology Core Facility at IPMB **.** 590 _Formation of Protein-DNA Complexes:_ _In vitro_ –translated proteins were immobilized on 591 HaloTag magnetic beads (Promega) by one hour incubation and washed three times with 592 PBS containing 0.005% NP-40 to remove unbound proteins. The purified TF–bead 593 complexes were then incubated with 100 ng of adapter-ligated genomic DNA and 10 μg 594 of Sheared Salmon Sperm DNA (Thermo Fisher) in binding buffer at room temperature 595 for 1 hour to allow specific protein–DNA interactions. Following incubation, beads were 596 magnetically separated and washed three times to eliminate nonspecific DNA. Bound 597 DNA was eluted by heating at 98°C for 10 minutes. 598 _DAP-seq analysis:_ Raw sequencing reads were trimmed using Trimmomatic (v0.39) [51] to 599 remove adapter sequences and low-quality bases, specifically employing a 600 SLIDINGWINDOW:4:15 quality filter and retaining reads with a minimum length of 36 bp. 601 Trimmed reads were aligned to the _Marchantia_ reference genome (v6.1r2) or _Arabidopsis_ 602 reference genome (TAIR 10) using Bowtie2 (v2.4.4) [52] in end-to-end mode with default 603 parameters. Significant binding regions were identified using MACS2 (v2.2.7.1) [53] against 604 a corresponding control sample, employing default settings. Peaks were retained if the 605 false discovery rate (FDR) with Benjamini-Hochberg correction was < 0.05. MpHSFB1 or 606 AtHSFB4 peak summits were assigned to target genes using the criterion as follows: a 607 peak was associated with a gene if its summit overlapped the proximal promoter region 608 (-1,000 bp to +500 bp relative to the annotated TSS). A Halo-tagged GFP was used for 609 negative control. For motif analysis, a 100 bp DNA sequence centered on each peak 610 summit (±50 bp) was extracted using BEDTools getfasta [54] . These sequences were 611 analyzed for enriched motifs using MEME [55], allowing for motif widths ranging from 6 to 612 20 base pairs with default background correction. 613 **Luciferase assay** 21 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 614 _Plasmid construction:_ Dual-luciferase reporter constructs were assembled via Golden 615 Gate cloning, with each binary vector containing three transcriptional units at defined 616 positions. Specifically, position P1 carried either a CaMV 35S promoter-driven MpHSFB1 617 6×HA effector cassette or an empty vector control; position P2 contained a firefly 618 luciferase (LUC2) reporter cassette driven by the candidate bHLH promoter (or a 35S 619 minimal promoter fused to the bHLH promoter); and position P3 contained a CaMV 35S 620 promoter-driven Renilla luciferase cassette to serve as an internal normalization control. 621 All bHLH promoters were amplified from Tak-1 genomic DNA using the specific primers 622 listed in Supporting Dataset S11. 623 _Agrobacterium-mediated transient expression in Nicotiana benthamiana:_ Binary vectors 624 were electroporated into _Agrobacterium tumefaciens_ strain GV3101. Individual colonies 625 were cultured for 2 days in selective LB medium at 28°C with shaking, pelleted by 626 centrifugation, and resuspended in infiltration buffer (10 mM MgCl₂, 10 mM MES, pH 5.6, 627 and 200 μM acetosyringone) to an OD₆₀₀ of 0.1. Suspensions were incubated at room 628 temperature for at least 30 minutes before being infiltrated into the abaxial side of 4-week 629 old _Nicotiana benthamiana_ leaves using a needleless syringe. For each construct, 630 infiltrations were performed on at least three independent leaves across three separate 631 plants. At 48 hours post-infiltration, leaf tissue was harvested for luciferase measurements. 632 _Dual-Glo luciferase assay:_ Infiltrated _N. benthamiana_ leaf tissue was collected from each 633 infiltration area, immediately frozen in liquid nitrogen, and homogenized in 150 μL Passive 634 Lysis Buffer (E1941, Promega). Lysates were clarified by centrifugation at 12,000 _g_ 1 5– 635 10 min at 4°C, and the supernatants were transferred to a white 96-well plate.Firefly and 636 Renilla luciferase activities were measured using the Dual-Glo® Luciferase Assay System 637 (E2920, Promega) according to the manufacturer’s instructions. Briefly, 25 μL of clarified 638 leaf extract was mixed with an equal volume of Dual-Glo® Luciferase Reagent. After 639 incubation for 10 min at room temperature, firefly luminescence was measured using a 640 plate reader (SpectraMax i3X, Molecular Devices). Subsequently, an equal volume of 641 Dual-Glo® Stop & Glo® Reagent was added to each well. After incubation for 10 min at 642 room temperature, Renilla luminescence was measured. Two technical replicates were 643 averaged before statistical analysis. Each independently infiltrated leaf area was treated 644 as one biological replicate. 22 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 645 **Generation of single-nuclei RNA-seq (snRNA-seq)** 646 _Nuclei isolation and snRNA-seq procedure:_ A total of 2.0 grams of 7-day-old _Marchantia_ 647 plants were chopped gently using an ice-cold razor blade in 2 mL of ice-cold modified 648 Galbraith buffer (45 mM MgCl2, 20 mM MOPS, 30 mM sodium citrate, 0.5% Triton X-100, 649 0.2 U/mL RNaseOUT™ Recombinant Ribonuclease Inhibitor, and 2 mM/mL DTT, pH 7.0) 650 for 4 minutes in a chilled petri dish, followed by incubation in a 4°C fridge for 10 minutes. 651 The lysate was subsequently filtered sequentially through 70 µm and 40 µm cell strainers 652 and transferred into an ice-cold 50 mL Falcon tube. Nuclei were pelleted at 500 × _g_ at 653 4°C for 5 minutes, after which the supernatant was carefully removed. The nuclei pellet 654 was gently resuspended in 1 mL of modified LB01 buffer (15 mM Tris pH 7.5, 0.5 mM 655 Spermine, 10 mM KCl, 20 mM NaCl, 0.1 U/mL RNaseOUT™ Recombinant Ribonuclease 656 Inhibitor, and 2 mM/mL DTT) and filter through 30 µm cell strainer into a 5 mL (12 × 75 657 mm) round-bottom tube followed by 2 µg/ml DAPI staining. 658 Flow cytometry was performed using a Beckman CytoFLEX SRT equipped with a 100 659 µm nozzle and operated at a low flow rate (psi < 20). The photomultiplier settings were 660 as follows: Forward Scatter (FSC) at 25 V, Side Scatter (SSC) at 25 V, and V450 (DAPI 661 channel) at 100 V. Sorting gates were empirically defined using FSC and SSC to exclude 662 debris and potential doublets. To establish baseline fluorescence, FACS analysis was 663 conducted on unstained control samples before sorting DAPI-stained nuclei. For nuclei 664 isolation, sort gates were set to specifically collect DAPI-positive nuclei. A total of 400,000 665 sorted nuclei were collected into 200 µL of 1X PBS + 1% BSA (0.4 U/mL RNaseOUT™ 666 Recombinant Ribonuclease Inhibitor) in 1.5 mL Eppendorf tubes within 30 minutes. The 667 nuclei were then pelleted by centrifugation at 500 × _g_ for 10 minutes at 4°C. The 668 supernatant was carefully removed, leaving approximately 80 µL. The quality and yield of 669 nuclei were checked under EVOS M7000 with C-Chip using the DAPI channel. About 670 16,000-19,000 nuclei of each sample were loaded into the 10X Genomics chip for cell 671 barcoding. The libraries were prepared following the manufacturer’s instructions with a 672 10X Chromium Single Cell 3’ Solution v3.1 kit. The paired-end libraries were sequenced 673 with the Illumina NextSeq2000 instrument at the NGS Core of BRCAS at Academia Sinica. 674 **snRNA-seq data analysis** 23 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 675 _Data processing and integration:_ snRNA-seq data from _M. polymorpha_ (Tak-1 plants and 676 Mp _hsfb1_ mutant) were processed using Seurat (v5.0.0) [56,57] . Raw gene-barcode matrices 677 generated by the 10X Genomics Cell Ranger (v9.0.0) pipeline were imported into R using 678 the Read10X function and converted to Seurat objects. Low-quality nuclei were filtered 679 based on the following criteria: ≥ 250 detected genes, ≥ 500 total unique molecular 680 identifiers (UMIs), and a log 10 ratio of genes per UMI > 0.8. Following filtering, each 681 dataset was normalized (NormalizeData), highly variable genes were identified 682 (FindVariableFeatures), and expression values were scaled (ScaleData). Principal 683 component analysis (PCA) was performed, and significant components were selected 684 based on elbow plots. Nuclei were clustered using a shared nearest-neighbor graph 685 (FindNeighbors, FindClusters) and visualized via UMAP (RunUMAP). 686 Potential multiplets were removed using DoubletFinder (v3.0) with an estimated 7.5% 687 doublet rate [58] . The optimal pK value was selected from the maximum mean-BCmetric 688 score across parameter sweeps, and adjustments were made for homotypic doublets 689 based on cluster composition. Only singlets were retained for downstream analysis. 690 Biological replicates for each genotype (Tak-1 and Mp _hsfb1_ ) were merged using merge 691 and flattened with JoinLayers to create unified RNA assays. To minimize batch effects, 692 Harmony integration [59] (RunHarmony) was performed on the first 30 PCs. These 693 integrated embeddings were used for subsequent clustering and visualization. 694 Finally, the integrated Tak-1 and Mp _hsfb1_ objects were aligned. Homologous nuclei were 695 mapped using FindTransferAnchors and TransferData, and cluster correspondences in 696 Mp _hsfb1_ were validated using SingleR with the Tak-1 dataset as a reference [60] . Nuclei 697 with a prediction confidence over 0.8 were retained. The final merged object, comprising 698 18,073 Tak-1 and 14,857 Mp _hsfb1_ nuclei, was batch-corrected again using Harmony 699 (grouping by genotype) and visualized by UMAP. 700 _Marker gene identification and functional annotation:_ Cluster-specific marker genes were 701 identified in the integrated Tak-1 and Mp _hsfb1_ dataset using FindAllMarkers (Wilcoxon 702 rank-sum test; min.pct = 0.25, logfc.threshold = 0.25, p_adj < 0.01). For each gene, 703 expression specificity was computed as the difference in detection frequency between 704 the target cluster and all other clusters. Marker lists were cross-referenced with a curated 24 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 705 _Marchantia_ genome (v6.1r2) to annotate transcription factors (TFs). To assign biological 706 functions, cluster markers were converted to _Arabidopsis_ orthologs via reciprocal BLAST 707 mapping. GO enrichment analysis was then performed using clusterProfiler (v4.10) with 708 _A. thaliana_ TERM2GENE and TERM2NAME mappings [61] . Significantly enriched 709 biological-process terms (FDR < 0.05) were visualized as heatmaps showing GO term 710 count and significance (−log 10 FDR). Cluster annotation was further refined by matching 711 marker genes to a previously published _Marchantia_ single-cell reference [26] . Marker 712 overlap ratios (matched / total annotated markers) were summarized, and the final 713 curated cluster identities were transferred to the Seurat object metadata. 714 _Analysis of differentially expressed genes:_ After integration of Tak-1 and Mp _hsfb1_ nuclei, 715 UMAP visualizations and cluster bar plots were generated to assess genotype-specific 716 cell distributions. Cell-type proportions between genotypes were quantified and compared 717 using speckle (v1.1.3) with a logit transformation (propeller) [62] . Differential expression 718 analysis between genotypes was conducted for each annotated cell type using 719 FindMarkers (Wilcoxon test). Genes with an adjusted p < 0.05 and |log 2 fold-change| > 720 0.25 were considered significant. DE results were summarized by gene counts and 721 proportions per cell type. To compare transcriptional programs, Jaccard similarity 722 coefficients were calculated between the Tak-1 and Mp _hsfb1_ cluster marker sets and 723 visualized as distance heatmaps. Finally, TF expression patterns were summarized by 724 their average log 2 fold-change across matched cell types. 725 _Trajectory inference:_ Pseudotime trajectories were reconstructed using Monocle3 726 (v1.3.1) [63-65] and Slingshot (v2.8.0) [66] on the integrated Seurat object to capture 727 developmental progression. UMAP coordinates served as the input for trajectory learning; 728 Monocle3 was used to identify the principal graph structure, while Slingshot refined 729 lineage relationships and root assignments. To identify genes dynamically associated with 730 pseudotime, generalized additive models (GAMs) were fitted using tradeSeq (v1.12.0) [67] . 731 The fitGAM function was applied for each lineage from Slingshot, and significant genes 732 were identified via association tests (p < 0.05). Smoothed expression profiles 733 (predictSmooth) were also generated at 100 evenly spaced pseudotime points for each 734 lineage. Genes consistently identified as pseudotime-dependent by both methods were 735 considered trajectory-relevant and used for regulatory network inference. 25 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 736 _Gene regulatory network (GRN) inference:_ To ensure regulatory network reconstruction 737 incorporated dynamic genes, the list of pseudotime-associated genes (from Monocle3 + 738 tradeSeq) was merged with expression-based filters before running SCENIC [68] . 739 Specifically, genes expressed in ≥ 1% of nuclei and with total expression > 3% of the 740 mean nucleus count were retained. The union of these filtered genes was the input 741 expression matrix for SCENIC. Regulatory networks (regulons) were inferred 742 independently for Tak-1 and Mp _hsfb1_ using SCENIC (v1.3.1) with motif enrichment 743 performed via RcisTarget (v1.29.0) [68] . Gene correlations were computed (runCorrelation), 744 and regulatory links were inferred using GENIE3 (runGenie3(exprMat_filtered, nParts = 745 20)), partitioning the dataset to optimize memory. This integrated approach ensured that 746 the inferred SCENIC networks reflected genes dynamically modulated along 747 developmental trajectories. 748 Regulon activity (AUC) matrices in Tak-1 or Mp _hsfb1_ SCENIC network were extracted 749 and aligned with cell-type annotations. Average regulon activities and regulon-specificity 750 scores (RSS) per cell type were computed. Genotype-specific and shared regulons were 751 compared using paired Wilcoxon tests, and differential activity was summarized in ΔZ 752 matrices (Tak-1 – Mp _hsfb1_ ). Network topology metrics—degree, betweenness, closeness, 753 eigenvector, hub, PageRank, and Katz centralities—were computed using igraph (v2.0.3). 754 Global comparisons were performed by Z-score normalization across genotypes, and 755 metric differences were visualized as dot heatmaps, where color represented the ΔZ 756 value and circle size corresponded to magnitude rank. All plots were generated using 757 ggplot2 (v3.5.1) or ComplexHeatmap (v3.21.0) [69], and all statistical analyses were 758 performed in R (v4.3.1). 759 **Generation of bulk RNA-seq and bulk RNA-seq analysis** 760 _Sample collection:_ For each sample, approximately 200 mg of gemmalings (from 3-, 7-, 761 and 14-day-old plants) or gemmalings with or without prolonged HS treatment (30 °C, 14 762 days) were promptly transferred to 2 ml Eppendorf tubes, flash-frozen in liquid nitrogen, 763 and stored for subsequent RNA extraction. Total RNA was extracted from the gemmalings 764 using the RNeasy Plant Mini Kit (Qiagen, Hilden, Germany). Library preparation, including 765 polyA enrichment, and sequencing of 150-bp paired-end reads were carried out on the 26 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 766 NovaSeq platform at the Genomic Technology Core (IPMB, Academia Sinica, Taiwan), 767 and GENOMICS (Taipei, Taiwan). 768 _Data analysis:_ Raw sequencing reads from _Marchantia_ were processed to remove 769 adaptors using fastp (v0.24.0) [70] . The cleaned reads were then mapped to the reference 770 genome ( _Marchantia_ genome v6.1r2) with STAR (v2.7.9a) using default parameters. 771 Differentially expressed genes (DEGs) were identified using the R package DESeq2 [71], 772 with pairwise comparisons performed against the corresponding genotype at day 3. In the 773 case of prolonged HS, the pairwise comparisons were performed against the 774 corresponding genotype at the control condition. DEGs were defined by a |Log 2 -fold 775 change| ≥ 1.0 and a False Discovery Rate (FDR) < 0.05, using the “lfcShrink ashr” method. 776 Z-scores of expressions were calculated from VST-transformed values. Time-course 777 differential expression analysis was conducted using the splineTimeR (v1.58.0) [72] 778 package to identify genes with dynamic expression changes across developmental 779 stages. DEGs were identified using the splineDiffExprs function with a smoothing degree 780 of freedom (df) = 2 and a significance threshold of FDR < 0.01. Separate pairwise 781 analyses were performed for Mp _hsfb1_ versus Tak-1 and Mp _HSFBox_ versus Tak-1 to 782 determine genotype-specific temporal expression patterns. Genes showing significant 783 temporal differences were exported as DEG lists for downstream analyses. Gene 784 Ontology (GO) enrichment analysis was performed using clusterProfiler (v4.10) [61] . 785 Significantly enriched GO terms (FDR < 0.05) were summarized using REVIGO [73] to 786 reduce redundancy. 787 **Cell-type-resolved co-expression network analysis using hdWGCNA** 788 _Data sources and preprocessing:_ Publicly available single-cell and single-nucleus 789 transcriptomic datasets from Arabidopsis, tomato, rice, and maize (Supporting Dataset 790 S8) were analyzed in R using Seurat for data handling and hdWGCNA (v0.4.0) [74] for cell 791 type-resolved co-expression network analysis. Genes were selected for network 792 construction based on expression prevalence to reduce sparsity and technical noise. To 793 improve the robustness of correlation estimates in single-cell data, metacells were 794 generated within each cell type and biological replicate by aggregating transcriptionally 27 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 795 similar cells in a low-dimensional space using k = 25 nearest neighbors. Normalized 796 metacell expression profiles were used as input for network inference. 797 _Cell-type-specific co-expression network construction:_ Co-expression networks were 798 constructed independently for each cell type using a signed weighted network framework 799 implemented in hdWGCNA. Soft-thresholding powers were selected based on scale-free 800 topology criteria, requiring a model fit of R² > 0.8. Adjacency matrices were transformed 801 into topological overlap matrices, followed by hierarchical clustering to identify co 802 expressed gene modules. Module eigengenes were calculated to summarize module 803 level expression patterns and harmonized across biological replicates where applicable. 804 _Quantification of HSFB connectivity:_ Module connectivity was quantified using 805 eigengene-based connectivity (kME) and node degree for each gene within its assigned 806 module. Genes assigned to the unclassified grey module were excluded from 807 downstream analyses. HSFB-associated modules were defined as modules containing 808 at least one HSFB gene within a given cell type and tissue and were interpreted as 809 putative HSFB co-expressed gene sets. To evaluate the role of HSFBs in each cell type, 810 their connectivity ranks were calculated relative to other transcription factors within the 811 same module and cell type based on kME values. Functional characterization of HSFB 812 associated modules was further supported by Gene Ontology enrichment analysis using 813 the full hdWGCNA gene set as the background. 814 **Integration of DAP-seq and scRNA-seq or snRNA-seq datasets** 815 _Acquisition and preprocessing of DAP-seq datasets:_ TF binding profiles were obtained 816 from two major DAP-seq resources: GSE60141 [33] and GSE299028 [32] . DAP-seq peaks 817 were assigned to target genes if the summit was located within the proximal promoter 818 region, consistently defined as -1000 bp upstream to +500 bp downstream of the 819 annotated transcription start site (TSS). The resulting TF-target lists were harmonized 820 and filtered for unique _Arabidopsis_ gene identifiers. 821 _Acquisition and preprocessing of scRNA-seq datasets:_ Single-cell and single-nucleus 822 transcriptomic data were sourced from four _Arabidopsis_ atlases: the root atlas [75] 823 (GSE152766) and the seedling, leaf, and flower bud atlases [32] (GSE297576). Seurat 824 objects were restricted to the RNA assay and normalized using the default log 28 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 825 normalization method (NormalizeData, LogNormalize). Cell-type annotations were 826 adopted from the original publications. 827 _Quantification_ _of_ _TF_ _activity_ _(AUCell):_ TF activity was quantified using 828 the AUCell framework [68], which ranks all genes by normalized expression per cell and 829 computes the area under the recovery curve (AUC) for each DAP-seq gene set. The AUC 830 reflects the enrichment of a TF’s target genes among the most highly expressed genes 831 in that cell. AUC scores were computed from the normalized data layer of each Seurat 832 object using default parameters in AUCell_buildRankings() and AUCell_calcAUC(), 833 producing a cell × TF activity matrix. 834 _Cell-type-specific enrichment analysis:_ TF activity enrichment was identified by 835 comparing AUC scores between cells of a given type and all other cells. Cohen’s _d_ (with 836 pooled standard deviation) was calculated using the _effectsize_ package to estimate effect 837 size, and one-sided Wilcoxon rank-sum tests assessed directional significance 838 (https://github.com/easystats/effectsize). P-values were adjusted using the Benjamini– 839 Hochberg procedure, and TFs with _d_ - 0 and FDR < 0.05 were considered significantly 840 enriched in the corresponding cell type. Results were visualized using dot 841 heatmaps generated 842 **Data availability** 843 Raw sequencing data have been deposited into the National Center for Biotechnology 844 Information Gene Expression Omnibus and will be made open upon publication. 845 **Acknowledgment** 846 We thank Ms. Shu-Jen Chou and Ms. Mei-Jen Fang from the Genomics Core at IPMB for 847 preparing the DAP-seq libraries and performing Sanger sequencing, respectively. We are 848 grateful to Dr. Meiyeh Jade Lu (NGS Core Facility, BRC, AS) for generating the snRNA 849 seq libraries, Dr. Wann-Neng Jane (Cell Biology Core, IPMB) for performing the SEM 850 experiments, and Dr. Yu-Chih Lin (Metabolomics Core, ABRC, AS) for measuring 851 cytokinin concentrations. We also acknowledge the Flow Cytometry Cores at IPMB and 852 IBMS, AS, for assistance with nuclei sorting, as well as the Bioinformatics Core at IPMB 853 for high-performance computing support. Furthermore, we thank Drs. Chih-Hang Wu and 29 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 854 Shu-Hsing Wu (IPMB, AS) for providing plasmids for the PPI and LUC reporter assay, 855 and Dr. Yucheng Long (National University of Singapore) for helpful discussions on 856 quantifying solidity using ImageJ. 857 **Funding** 858 This study is supported by funding from the National Science and Technology Council 859 (NSTC) 115-2311-B-001-027-MY3 and Academia Sinica (AS) AS-CDA-113-L03 provided 860 to TYW. 861 **Author’s contribution** 862 TYW conceptualized and oversaw all the experiments and analyses. HCT and KHL 863 designed and performed all the experiments with the assistance of HJY and TYW. YRL 864 and MJL performed PPI experiments. LT performed snRNA-seq and downstream analysis. 865 MRY and TYW performed formal bioinformatic analyses with assistance from CYC. TYW 866 wrote the manuscript with assistance from all authors. HCT, KHL, CYC, and TYW edited 867 the manuscript. 868 **Competing interests** 869 The authors declare no competing interests. 870 **Figure Legend** 871 **Figure 1. MpHSFB1 is involved in terrestrial adaptation.** **A.** Schematic illustration of 872 four key features associated with adaptation to terrestrial life. 1. Plants evolved rhizoid 873 structures for water and nutrient absorption. 2. They also developed air pore structures 874 for gas exchange. 3–4. In addition, hormonal and ROS signaling pathways evolved to 875 enable responses to unfavorable environments. B. Representative images of Tak-1, 876 Mp _hsfb1_, and Mp _HSFB1ox_ plants from day 0 to day 30. **C.** Representative images 877 showing developmental traits at day 14. **D.** Heatmaps showing quantitative analyses of 878 developmental traits. Z-scores were calculated across developmental stages and 879 genotypes. p-values were determined using Student’s t-test between mutants and the 880 wild type. Two independent experiments were performed with _n_ = 30–60 per time point. 881 The raw measurements and complete statistical results are provided in Supporting 882 Dataset S1. **E.** Schematic overview of the datasets collected in this study. 30 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 883 **Figure 2. snRNA-seq reveals cell-type changes in Mp** _**hsfb1**_ **mutants. A.** UMAP 884 visualization of major cell-type clusters in Tak-1 plants and Mp _hsfb1_ mutants. **B.** Bar graph 885 showing cell-type enrichments between Tak-1 and Mp _hsfb1_ . **p < 0.01, hypergeometric 886 test. **C.** Heatmap showing cell-type similarity between Tak-1 and Mp _hsfb1_ . The color key 887 represents Jaccard distance, where lower values indicate higher similarity. **D.** Ridgeline 888 plot showing pseudotime distributions across cell types in Tak-1 and Mp _hsfb1_ . **p < 0.01, 889 empirical cumulative distribution function (ECDF), Kolmogorov–Smirnov test. **E.** Heatmap 890 showing transcription factor (TF) families differentially expressed (DE) between Tak-1 and 891 Mp _hsfb1_, enriched in high (>10%) or low (<10%) proportions of cells. Percentages were 892 calculated as the number of DE TFs over the total TFs within each family. Note that only 893 the TF families with a number of DE TFs greater than 2 were shown. The full list was 894 presented in the Supporting Dataset S2. 895 **Figure 3. The Mp** _**hsfb1**_ **regulons underwent functional rewiring through altered** 896 **gene targeting.** **A.** Comparison of regulon properties between Tak-1 and Mp _hsfb1_ . The 897 Venn diagram shows the number of overlapping regulons between the two genotypes. **B.** 898 Table comparing the network composition of Tak-1 and Mp _hsfb1_ . **C.** Dot plot showing the 899 enrichment (hypergeometric test) of TF centrality metrics between Tak-1 and Mp _hsfb1_ . 900 **D.** Dot plots represent the scaled regulon specificity score (RSS) after z-score 901 transformation in Tak-1 and Mp _hsfb1_ . Dot size indicates the Δz-score between Tak-1 and 902 Mp _hsfb1_, and dot color denotes the cell type with the maximum RSS for each regulon. **E.** 903 Network visualization showing representative GO terms from each regulon in selected 904 cell types. Edge width represents the significance of GO term enrichment, and node size 905 indicates the number of connected regulons (TFs). 906 **Figure 4.** **MpHSFB1 regulates cytokinin homeostasis during the early** 907 **developmental stage of** _**Marchantia**_ **. A.** Hierarchical clustering of the time-course bulk 908 RNA-seq dataset. The line represents the mean gene expression of each genotype group. 909 **B.** Heatmap showing the top three significantly enriched GO (Biological Process, BP) 910 terms from each cluster (q-value < 0.01). **C.** Line graphs showing cis-zeatin levels. Two 911 biological replicates with two technical replicates each were performed for each genotype 912 and condition. *p < 0.05 was determined by a pairwise Student’s t-test. **D.** Representative 913 photos of cytokinin treatments in Tak-1, Mp _hsfb1_, and Mp _HSFB1ox_ plants. Note that the 31 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 914 photo for Tak-1 plants is from the same set as in Figure 6A (control condition), as the 915 experiments were performed simultaneously. **E.** Boxplots showing thallus area for each 916 genotype and condition. Dots represent individual plant data; p-values were determined 917 by the Mann-Whitney U test with Bonferroni correction ( _n_ = 60 plants per line). Boxplots 918 display all data points (dots), median values, first and third quartiles, and whiskers 919 indicating minimum and maximum values. Fold change was calculated as the ratio of the 920 mean thallus area between control and cytokinin treatments for each genotype. **F.** 921 Heatmap showing the enrichment (hypergeometric distribution) of TF families in each 922 cluster. Only TF families enriched in at least one cluster are shown. **G.** Integrated 923 regulatory network of Mp _HSFB1_ -mediated Mp _bHLH_ regulation. The network integrates 924 data from time-course RNA-seq clusters, public RNA-seq datasets, snRNA-seq datasets, 925 regulons, DAP-seq, and known mutant phenotypes. Node size corresponds to the 926 number of lines of supporting evidence (bulk RNA-seq, snRNA-seq, and DAP-seq). Solid 927 lines represent physical DNA–protein or protein–protein binding interactions. Edge colors 928 indicate whether the regulation is associated with a rhizoid- or gemma-phenotype-related 929 bHLH gene. See Supporting Datasets S6 and S7 for detailed information. 930 **Figure 5.** **MpbHLHs regulate cytokinin homeostasis during the early developmental** 931 **stage of** _**Marchantia**_ **. A.** Representative images of Mp _bhlh_ knockout mutants at day 14. 932 **B.** Boxplot showing the rhizoid-to-thallus ratio. P-values were determined by one-way 933 ANOVA followed by Tukey HSD post hoc test ( _n_ = 10-12 plants per line). **C.** Boxplot 934 showing the _cis_ -zeatin levels. P-values were determined by one-way ANOVA followed by 935 Tukey HSD post hoc test ( _n_ = 6-8 independent measurements). **D.** Scatter plot showing 936 the Spearman correlation between the rhizoid/thallus ratio and cytokinin concentration. 937 Each data point represents the mean value, with error bars indicating the standard error 938 (S.E.). **E.** Representative photos of cytokinin treatments in Tak-1, Mpb _hlh21-14_, 939 Mp _bhlh21-16_, Mp _bhlh40-5_, and Mp _bhlh40-6_ plants. **F.** Boxplots showing thallus area for 940 each genotype and condition. P-values were determined by the Mann-Whitney U test with 941 Bonferroni correction ( _n_ = 30-60 plants per line). Fold change was calculated as the ratio 942 of the mean thallus area between control and cytokinin treatments for each genotype. **G.** 943 **& H.** UpSet plots and heatmaps showing genotype-specific DEGs of Mp _bhlh21_ ( **G** ) and 944 Mp _bhlh40_ ( **H** ) mutants at days 7 and 14 relative to wild-type Tak-1 plants. For all panels, 32 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 945 boxplots display individual data points (dots), median values, first and third quartiles, and 946 whiskers indicating minimum and maximum values. 947 **Figure 6.** **MpHSFB1 regulates prolonged HS tolerance via the ROS response. A.** 948 Representative top-view and side-view images after prolonged HS treatment at day 14. 949 **B.** Heatmaps showing quantitative analyses of developmental traits. Z-scores were 950 calculated across developmental stages and genotypes. p-values were determined using 951 Student’s t-test between mutants and the wild type. Two independent experiments were 952 performed with _n_ = 30–60 per time point. Raw measurements and complete statistical 953 results are provided in Supporting Dataset S1. **C.** Boxplots showing DAB staining 954 intensities in each genotype and condition. The y-axis denotes grayscale pixel values (0 955 [black]–255 [white]). P-values were determined by the Mann-Whitney U test with 956 Bonferroni correction ( _n_ = 12–60 plants per line). **D.** UpSet plot showing up- and down 957 regulated DEGs under prolonged HS. The red bars denote Mp _hsfb1_ - or Mp _HSFB1ox_ 958 specific genes subjected to further analysis. **E.** Heatmap showing representative enriched 959 GO (BP) terms (FDR < 0.05) from the groups in (D). **F.** Heatmap showing log₂FC values 960 of selected TFs (|∆( _B1ox_ - _b1_ )| ≥ 1) across genotypes under prolonged HS, highlighting 961 expression differences between Mp _HSFB1ox_ lines and Mp _hsfb1_ mutants. Asterisks 962 denote DE TFs, and red labels indicate TFs selected for further validation. **G.** 963 Representative images of Mp _bhlh51_ and Mp _nac6_ knockout mutants at day 14 under 964 prolonged HS. **H.** Boxplots showing thallus area for each genotype and condition. P 965 values were determined by the Mann-Whitney U test with Bonferroni correction ( _n_ = 60 966 plants per line). **I.** Dot plot showing the grayscale ratio under prolonged HS relative to 967 control conditions for each genotype. Error bars indicate the standard error (S.E.). A lower 968 ratio indicates a higher ROS level. *p < 0.05 was determined by a pairwise Student’s t 969 test. **J.** Venn diagram showing the number of genotype-specific DEGs in _Mpbhlh51_ and 970 _Mpnac6_ mutants under prolonged HS conditions. **K.** Significantly enriched GO BP terms 971 corresponding to the DEGs identified in (J). For all panels, boxplots display individual data 972 points (dots), median values, first and third quartiles, and whiskers indicating minimum 973 and maximum values. 974 **Figure 7.** **MpHSFB1 and its orthologs coordinate tissue- and cell type-specific** 975 **transcriptional programs across land plants. A.** Schematic illustrating the 33 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 976 experimental design for swapDAP-seq. **B & C.** Venn diagrams showing the overlapping 977 orthogroups identified by swapDAP-seq using MpHSFB1 ( **B** ) or AtHSFB4 ( **C** ) bound to 978 the Marchantia and Arabidopsis genomes. **D.** Heatmap showing the enriched GO terms 979 identified from genes bound exclusively in Arabidopsis, exclusively in Marchantia, or in 980 both species from ( **B** .). **E.** Heatmap showing the similarity of TF regulation across four 981 Arabidopsis tissues, calculated using the Pearson correlation coefficient. **F.** Proposed role 982 of MpHSFB1 in land plant adaptation. The model illustrates how MpHSFB1 coordinates 983 distinct transcriptional programs for development and stress resilience. In _Marchantia_, 984 MpHSFB1 cooperates with specific bHLH partners to regulate cytokinin homeostasis 985 during early development and ROS responses under prolonged heat stress. This 986 regulatory flexibility highlights the co-option of HSFB1 networks in orchestrating tissue 987 specific adaptations across land plants. 988 **Supplemental information** 989 **Supporting Figure S1. Characterization of Mp** _**hsfb1**_ **knockout and Mp** _**HSFB1ox**_ 990 **overexpression lines.** **A.** Bar graph showing MpHSFB1 expression levels in Tak-1 plants, 991 as well as Mp _HSFB1ox-1_ and Mp _HSFB1ox-6_ lines. Fold change was calculated relative 992 to endogenous MpHSFB1 expression in Tak-1. _N_ = 3 biological replicates. **B.** 993 Representative photos of the second independent lines of Mp _hsfb1_ and Mp _HSFB1ox_ 994 plants. Note that the photo for Tak-1 plants is from the same set as in Figure 1, as the 995 experiments were performed simultaneously. **C.** Phylogenetic tree of plant HSFB 996 sequences. The maximum-likelihood tree was constructed using the DNA-binding domain 997 (DBD) sequences of selected taxa. The HSFB1–B4 groups are color-coded, and all 998 sequences used for tree construction are listed in Supporting Dataset 1. 999 **Supporting Figure S2. Quality control of snRNA-seq data and identification of cell-** 1000 **type-specific marker genes.** **A.** Output files from the Cell Ranger software for each 1001 snRNA-seq sample. **B.** Dot plots showing the marker genes for each known cell type. 1002 Information regarding marker genes and cell types was taken from a previous study [26] . 1003 Dot size indicates the percentage of cells expressing the gene (% expressed). Dot colors 1004 represent the proportional amount of expression of each gene in each cell type, with 1005 deeper colors indicating higher expression levels. 34 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 1006 **Supporting Figure S3. Characterization of cell-type-specific DEGs in Tak-1 plants** 1007 **and Mp** _**hsfb1**_ **mutants. A.** Dot plot summarizing the proportion of DEGs found in each 1008 cell type in Tak-1 plants and Mp _hsfb1_ mutants. **B.** Dot plot showing the counts of 1009 differentially expressed transcription factors (DETFs) in each cell type. The color 1010 represents the odds ratio. **C.** Dot plot showing the overrepresented GO terms in each cell 1011 type. The size of the circle represents the DEG count. **D.** Representative images of oil 1012 body staining for each genotype. **E.** Boxplot showing the number of oil bodies per gemma 1013 for each genotype. Boxplots display all individual data points (dots), median values, first 1014 and third quartiles, and whiskers indicating minimum and maximum values. _P_ -values 1015 were determined by a pairwise Student’s _t_ -test. _n_ = 10. 1016 **Supporting Figure S4. Pseudotime trajectory and cell lineage analysis in snRNA-** 1017 **seq. A.** UMAP showing the lineage analysis from Slingshot. In total, seven possible 1018 lineages were identified. **B.** Heatmap showing the expression patterns of differentially 1019 expressed bHLH and NAC TFs in each lineage within the Tak-1 and Mp _hsfb1_ datasets. 1020 The color bar indicates the relative expression level. 1021 **Supporting Figure S5. snRNA-seq network analysis and identification of hormone-** 1022 **related gene expression.** **A.** Boxplots showing seven metrics of shared regulons 1023 identified from Tak-1 and Mp _hsfb1_ mutants. Dots represent individual plant data; p-values 1024 were determined by the Mann-Whitney U test with Bonferroni correction. Boxplots display 1025 all data points (dots), median values, first and third quartiles, and whiskers indicating 1026 minimum and maximum values. **B.** Dot plots showing the scaled regulon specificity score 1027 (RSS) after z-score transformation in Tak-1 and Mp _hsfb1_ for each network metric. Dot 1028 size indicates the value of Δz-score, and dot color denotes the Δz-score between Tak-1 1029 and Mp _hsfb1_ . **C.** Heatmap showing the genotype-specific scaled regulon specificity score 1030 (RSS) after z-score transformation in Tak-1 and Mp _hsfb1_ . The color bar represents the 1031 activity of regulons. **D & E.** Heatmaps showing relative expression of genes related to 1032 cytokinin ( **D** ) and abscisic acid (ABA) ( **E** ) biosynthesis and signaling in each annotated 1033 cell type and genotype. Gene expression levels are represented by aggregated unique 1034 molecular identifier (UMI) counts across all cells belonging to the same cell type within 1035 one genotype. The adjusted p-value was determined by the hypergeometric distribution. 35 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 1036 **Supporting Figure S6. Time-course bulk RNA-seq analysis of Mp** _**hsfb1**_ **,** 1037 **Mp** _**HSFB1ox**_ **, and Tak-1 plants.** **A.** Principal Component Analysis (PCA) comparing 1038 RNA-seq biological replicates from each genotype at days 3, 7, and 14. **B.** Bar graph 1039 displaying the number of differentially expressed genes (DEGs) in each genotype at days 1040 3, 7, and 14. Note that data from Tak-1 plants at day 3 were used as the baseline for 1041 pairwise comparisons between mutants and wild-type plants; therefore, no DEGs are 1042 shown for Tak-1 at day 3. For comparisons at day 7 and day 14, the respective genotype’s 1043 data at day 3 were used as the baseline. **C.** Heatmap showing the expression profiles of 1044 cytokinin biosynthesis genes in Tak-1, Mp _HSFB1-ox_, and Mp _hsfb1_ plants at days 7 and 1045 14 relative to day 3. Asterisks indicate significant DEGs (|log2FC| ≥ 1, FDR < 0.01). **D.** 1046 Top three binding motifs determined by DAP-seq and their closest known TFBSs. 1047 **Supporting** **Figure** **S7.** **Validation** **of** **MpHSFB1–MpbHLH** **protein-protein** 1048 **interactions in tobacco leaves.** **A.** Schematic illustrating construct designs for protein 1049 protein interaction using the split GAL4-RUBY assay and the arrangement of infiltration 1050 sites on tobacco leaves. The known interaction between bHLH16 and SMF served as the 1051 positive control [76] . **B.** Representative photos showing betalain production in tobacco 1052 leaves for each tested TF. Experiments were repeated three times with three leaves per 1053 replicate, yielding consistent results. **C. & D** . Independent biological replicates (second 1054 and third) of the experiments shown in A and B. Note that in the third biological replicate, 1055 leaves were destained with 95% ethanol to remove chlorophyll for better visualization. **E.** 1056 Dual-luciferase reporter activity (firefly/ _Renilla_ ) measured in _N. benthamiana_ co 1057 transfected with each _MpBHLHpro:firefly_ reporter and either vehicle (white bars) or the 1058 MpHSFB1 effector (grey bars). For each reporter, activity is shown under the reporter 1059 alone (Promoter only) and with the effector driven by the CaMV 35S promoter 1060 (35S:promoter). Bars represent mean ± SE (n = 4 biological replicates). Different letters 1061 denote statistically significant differences (p < 0.05, two-way ANOVA with Tukey's post 1062 hoc test). 1063 **Supporting Figure S8. Characterization of Mp** _**bhlh**_ **knockout mutants. A.** Indel 1064 sequences of Mp _bhlh21-14_, Mp _bhlh21-16_, Mp _bhlh28-18_, Mp _bhlh40-5_, Mp _bhlh40-6_, 1065 Mp _bhlh51-33_, Mp _bhlh51-96_, Mp _nac6-3_, and Mp _nac6-5_ plants. Red text indicates 1066 deletions, and yellow boxes highlight the gRNA target sites. **B.** Representative photos 36 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 1067 illustrating gemma cup morphology in Mp _bhlh28-18_, Mp _bhlh40-5_, and Mp _bhlh40-6_ plants 1068 compared to Tak-1 plants under control conditions. Photos were taken after 24 days of 1069 growth on ½ × B5 plates. 1070 **Supporting Figure S9. Characterization of Mp** _**bhlh**_ **/Mp** _**HSFB1ox**_ **lines.** **A.** 1071 Representative images of Mp _bhlh_ /Mp _HSFB1ox_ lines and their corresponding vector 1072 controls after 14 days of growth on ½ × B5 plates. **B.** Relative _MpHSFB1_ expression 1073 levels measured by RT-qPCR across all transgenic lines shown in ( **A** ). **C.–E.** Boxplots 1074 showing thallus area ( **C** ), solidity ( **D** ), and rhizoid-to-thallus ratio ( **E** ) for each genotype. 1075 Lowercase letters indicate statistically significant differences determined by one-way 1076 ANOVA with Tukey’s post hoc HSD test. For all panels, boxplots display individual data 1077 points (dots), median values, first and third quartiles, and whiskers indicating minimum 1078 and maximum values. 1079 **Supporting Figure S10. Analysis of bulk RNA-seq datasets in Mp** _**hsfb1**_ **,** 1080 **MpH** _**SFB1ox**_ **, and Tak-1 plants under prolonged HS. A.** Principal Component Analysis 1081 (PCA) comparing RNA-seq biological replicates from each genotype under each condition. 1082 **B.** Bar graph displaying the number of DEGs in each genotype under prolonged HS 1083 treatment. **C.** Representative photos of the second independent lines of Mp _hsfb1_ and 1084 Mp _HSFB1ox_ plants under prolonged HS treatment. Note that the photo for Tak-1 plants 1085 is from the same set as in Figure 5A (HS condition), as the experiments were performed 1086 simultaneously. **D.** Boxplots showing thallus area (in pixels) for each genotype and 1087 condition. P-values were determined by the Mann-Whitney U test with Bonferroni 1088 correction ( _n_ = 60 plants per line). Fold change was calculated as the ratio of the mean 1089 thallus area between control and HS treatments for each genotype. **E.** Dual-luciferase 1090 reporter activity (firefly/Renilla) measured in _N. benthamiana_ co-transfected with each 1091 _MpBHLHpro:firefly_ reporter and either vehicle (white bars) or the MpHSFB1 effector (grey 1092 bars). For each reporter, activity is shown under the reporter alone (Promoter only) and 1093 with the effector driven by the CaMV 35S promoter (35S:promoter). Bars represent mean 1094 ± SE (n = 4 biological replicates). Different letters denote statistically significant 1095 differences (p < 0.05, two-way ANOVA with Tukey's post-hoc test). **F.** Boxplots showing 1096 DAB staining intensities in each genotype and condition. The y-axis denotes grayscale 37 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 1097 pixel values (0 [black]–255 [white]). P-values were determined by the Mann-Whitney U 1098 test with Bonferroni correction ( _n_ = 60-100 plants per line). For all panels, boxplots display 1099 individual data points (dots), median values, first and third quartiles, and whiskers 1100 indicating minimum and maximum values. 1101 **Supporting Figure S11. Regulatory specificity of MpHSFB1-mediated MpbHLH** 1102 **regulation.** **A.** Boxplots showing thallus area for each genotype and condition. Dots 1103 represent individual plant data; _p_ -values were determined by the Mann-Whitney _U_ test 1104 with Bonferroni correction ( _n_ = 60 plants per line). **B.** (Left) Boxplots showing DAB staining 1105 intensities for each genotype and condition. The y-axis denotes grayscale pixel values 1106 from 0 (black) to 255 (white). P-values were determined by the Mann-Whitney _U_ test with 1107 Bonferroni correction ( _n_ = 12–60 plants per line). (Right) Dot plot showing the grayscale 1108 ratio under prolonged HS relative to control conditions for each genotype. Error bars 1109 indicate the standard error (S.E.), where a lower ratio indicates higher ROS levels. **C.** 1110 Boxplot showing the rhizoid-to-thallus ratio. P-values were determined by one-way 1111 ANOVA followed by Tukey's HSD _post hoc_ test ( _n_ = 10–12 plants per line). **D.** Boxplot 1112 showing _cis_ -zeatin levels. P-values were determined by one-way ANOVA followed by 1113 Tukey's HSD _post hoc_ test ( _n_ = 6–8 independent measurements). **E.** Boxplots showing 1114 thallus area for each genotype and condition. P-values were determined by the Mann 1115 Whitney _U_ test with Bonferroni correction ( _n_ = 30–60 plants per line). For all panels, 1116 boxplots display individual data points (dots), median values, first and third quartiles, and 1117 whiskers indicating minimum and maximum values. 1118 **Supporting Figure S12. HSFB co-expression networks identified by hdWGCNA** 1119 **across various land plants.** (A) Network showing 8 groups of GO term categories with 1120 their corresponding functions. The color of the node represents the significance of the GO 1121 biological process (BP) term (-log10(FDR)), and the size of the node represents the total 1122 number of annotated genes. (B) Bar graphs showing the proportion of GO terms in each 1123 GO category for selected HSFBs. P-values were determined by Fisher’s exact test 1124 (comparisons made against MpHSFB1). 1125 **Supporting Figure S13. Comparative swapDAP-seq analysis in** _**Marchantia**_ **and** 1126 _**Arabidopsis**_ **.** **A.** Metaplots showing the distribution of MpHSFB1 on _Marchantia_ and 38 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 1127 _Arabidopsis_ genomic DNA in the region (± 2 kb) surrounding transcription start sites 1128 (TSSs). The sequence logo of the binding motif identified from _Arabidopsis_ genomic DNA 1129 is shown in the bottom panel. **B.** Venn diagrams showing overlapping orthogroups 1130 identified by swap-DAP-seq using MpHSFA1 or AtHSFA1b bound to the _Marchantia_ and 1131 _Arabidopsis_ genomes. **C.** Dot plot showing overrepresented GO terms in each cell type. 1132 Circle size represents the count of target genes for each HSF. 1133 **Supporting Figure S14. Multi-tissue regulatory map of HSFB- and bHLH-mediated** 1134 **cell-type-specific gene expression programs.** **A–C.** Dot plots showing the enrichment 1135 or depletion of target gene expression for each HSFA and HSFB (horizontal axis) across 1136 labeled cell types in the Arabidopsis root ( **A** ), seedling ( **B** ), and flower bud ( **C** ) snRNA 1137 seq atlases (vertical axis). Dot size represents the FDR-corrected −log 10 (p-value) from a 1138 Wilcoxon rank-sum test of summarized target expression scores for each cell-type label 1139 versus all other cells. Dot color represents the effect size (Cohen’s d) of expression 1140 enrichment (red) or depletion (blue). Enrichments with FDR-corrected p < 0.01 are shown. 1141 **D.** Same as A–C but for the Marchantia snRNA-seq atlas. 1142 **Supporting Dataset S1. Raw phenotypic data under control, cytokinin, and** 1143 **prolonged HS conditions, and endogenous cytokinin measurements.** 1144 **Supporting Dataset S2. Quality control of snRNA-seq data, cell type annotation,** 1145 **and identification of differentially expressed genes.** 1146 **Supporting Dataset S3. Results of cell lineage analysis using Slingshot.** 1147 **Supporting Dataset S4. snGRN analysis in Tak-1 and Mp** _**hsfb1**_ **plants.** 1148 **Supporting Dataset S5. GO enrichment analysis results for each regulon.** 1149 **Supporting Dataset S6. List of DEGs from 3-, 7-, and 14-day-old Tak-1, Mp** _**hsfb1**_ **,** 1150 **and Mp** _**HSFB1ox**_ **plants, and list of MpHSFB1-targeted genes from DAP-seq.** 1151 **Supporting Dataset S7. Lists of DEGs from 7- and 14-day-old Mp** _**bhlh21**_ **and** 1152 **Mp** _**bhlh40**_ **mutants, alongside re-analyses of publicly available Mp** _**bHLH**_ **mutant** 1153 **RNA-seq datasets.** 1154 **Supporting Dataset S8. List of DEGs from Tak-1, Mp** _**hsfb1**_ **, Mp** _**HSFB1ox,**_ **Mp** _**bhlh51**_ 1155 _**and**_ **Mp** _**nac6**_ **plants under prolonged HS.** 39 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 1156 **Supporting Dataset S9. List of co-expression modules and corresponding GO** 1157 **terms identified by hdWGCNA in** _**Arabidopsis**_ **, rice, tomato, maize, and** _**Marchantia**_ 1158 **scRNA-seq or snRNA-seq datasets.** 1159 **Supporting Dataset S10. List of genes identified from swapDAP-seq and HSFB** 1160 **regulons identified from different scRNA-seq or snRNA-seq datasets.** 1161 **Supporting Dataset S11. List of primers and plasmids used in this study.** 1162 **References** 1163 1. de Vries, J., and Archibald, J.M. (2018). Plant evolution: landmarks on the path to 1164 terrestrial life. New Phytol _217_, 1428-1434. 10.1111/nph.14975. 1165 2. Furst-Jansen, J.M.R., de Vries, S., and de Vries, J. (2020). Evo-physio: on stress 1166 responses and the earliest land plants. J Exp Bot _71_, 3254-3269. 1167 10.1093/jxb/eraa007. 1168 3. Scheres, B., and van der Putten, W.H. (2017). The plant perceptron connects 1169 environment to development. Nature _543_, 337-345. 10.1038/nature22010. 1170 4. Powell, A.E., and Heyl, A. (2023). The origin and early evolution of cytokinin 1171 signaling. Front Plant Sci _14_, 1142748. 10.3389/fpls.2023.1142748. 1172 5. Kieber, J.J., and Schaller, G.E. (2018). Cytokinin signaling in plant development. 1173 Development _145_ . 10.1242/dev.149344. 1174 6. Hori, K., Maruyama, F., Fujisawa, T., Togashi, T., Yamamoto, N., Seo, M., Sato, S., 1175 Yamada, T., Mori, H., Tajima, N., et al. (2014). Klebsormidium flaccidum genome 1176 reveals primary factors for plant terrestrial adaptation. Nat Commun _5_, 3978. 1177 10.1038/ncomms4978. 1178 7. Rashotte, A.M. (2021). The evolution of cytokinin signaling and its role in 1179 development before Angiosperms. Semin Cell Dev Biol _109_, 31-38. 1180 10.1016/j.semcdb.2020.06.010. 1181 8. Li, S.-M., Zheng, H.-X., Zhang, X.-S., and Sui, N. (2021). Cytokinins as central 1182 regulators during plant growth and stress response. Plant Cell Reports _40_, 2711183 282. 10.1007/s00299-020-02612-1. 1184 9. Wu, T.Y., Hoh, K.L., Boonyaves, K., Krishnamoorthi, S., and Urano, D. (2022). 1185 Diversification of heat shock transcription factors expanded thermal stress 1186 responses during early plant evolution. Plant Cell _34_, 3557-3576. 1187 10.1093/plcell/koac204. 1188 10. Wang, X., Shi, X., Chen, S., Ma, C., and Xu, S. (2018). Evolutionary Origin, 1189 Gradual Accumulation and Functional Divergence of Heat Shock Factor Gene 1190 Family with Plant Evolution. Front Plant Sci _9_, 71. 10.3389/fpls.2018.00071. 1191 11. Zheng, L., Zhang, Q., Wang, C., Wang, Z., Gao, J., Zhang, R., Shi, Y., and Zheng, 1192 X. (2025). The heat shock factor HSFB1 coordinates plant growth and drought 1193 tolerance in Arabidopsis. The Plant Journal _121_, e17258. 1194 https://doi.org/10.1111/tpj.17258. 1195 12. Ikeda, M., Mitsuda, N., and Ohme-Takagi, M. (2011). Arabidopsis HsfB1 and 1196 HsfB2b act as repressors of the expression of heat-inducible Hsfs but positively 40 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 1197 regulate the acquired thermotolerance. Plant Physiol _157_, 1243-1254. 1198 10.1104/pp.111.179036. 1199 13. Pardal, R., Scheres, B., and Heidstra, R. (2023). SCHIZORIZA domain-function 1200 analysis identifies requirements for its specific role in cell fate segregation. Plant 1201 Physiol _193_, 1866-1879. 10.1093/plphys/kiad456. 1202 14. Aki, S.S., Nishihama, R., Kohchi, T., and Umeda, M. (2019). Cytokinin signaling 1203 coordinates development of diverse organs in Marchantia polymorpha. Plant 1204 Signal Behav _14_, 1668232. 10.1080/15592324.2019.1668232. 1205 15. Aki, S.S., Mikami, T., Naramoto, S., Nishihama, R., Ishizaki, K., Kojima, M., 1206 Takebayashi, Y., Sakakibara, H., Kyozuka, J., Kohchi, T., and Umeda, M. (2019). 1207 Cytokinin Signaling Is Essential for Organ Formation in Marchantia polymorpha. 1208 Plant Cell Physiol _60_, 1842-1854. 10.1093/pcp/pcz100. 1209 16. Fu, J., Zhou, C., Ma, F., Zhao, J., Yu, F., and Cui, H. (2024). The PLETHORA 1210 Homolog in Marchantia polymorpha is Essential for Meristem Maintenance, 1211 Developmental Progression, and Redox Homeostasis. Plant and Cell Physiology 1212 _65_, 1231-1244. 10.1093/pcp/pcae055. 1213 17. Gao, F., and Dubos, C. (2024). The arabidopsis bHLH transcription factor family. 1214 Trends in Plant Science _29_, 668-680. https://doi.org/10.1016/j.tplants.2023.11.022. 1215 18. Pires, N., and Dolan, L. (2009). Origin and Diversification of Basic-Helix-Loop-Helix 1216 Proteins in Plants. Mol Biol Evol _27_, 862-874. 10.1093/molbev/msp288. 1217 19. Catarino, B., Hetherington, A.J., Emms, D.M., Kelly, S., and Dolan, L. (2016). The 1218 Stepwise Increase in the Number of Transcription Factor Families in the 1219 Precambrian Predated the Diversification of Plants On Land. Mol Biol Evol _33_, 1220 2815-2819. 10.1093/molbev/msw155. 1221 20. Kanaoka, M.M., Pillitteri, L.J., Fujii, H., Yoshida, Y., Bogenschutz, N.L., 1222 Takabayashi, J., Zhu, J.-K., and Torii, K.U. (2008). SCREAM/ICE1 and SCREAM2 1223 Specify Three Cell-State Transitional Steps Leading to Arabidopsis Stomatal 1224 Differentiation The Plant Cell _20_, 1775-1785. 10.1105/tpc.108.060848. 1225 21. Bruex, A., Kainkaryam, R.M., Wieckowski, Y., Kang, Y.H., Bernhardt, C., Xia, Y., 1226 Zheng, X., Wang, J.Y., Lee, M.M., Benfey, P., et al. (2012). A gene regulatory 1227 network for root epidermis cell differentiation in Arabidopsis. PLoS Genet _8_, 1228 e1002446. 10.1371/journal.pgen.1002446. 1229 22. Karas, B., Amyot, L., Johansen, C., Sato, S., Tabata, S., Kawaguchi, M., and 1230 Szczyglowski, K. (2009). Conservation of Lotus and Arabidopsis Basic Helix-Loop1231 Helix Proteins Reveals New Players in Root Hair Development Plant Physiology 1232 _151_, 1175-1185. 10.1104/pp.109.143867. 1233 23. Proust, H., Honkanen, S., Jones, Victor A.S., Morieri, G., Prescott, H., Kelly, S., 1234 Ishizaki, K., Kohchi, T., and Dolan, L. (2016). RSL Class I Genes Controlled the 1235 Development of Epidermal Structures in the Common Ancestor of Land Plants. 1236 Curr Biol _26_, 93-99. https://doi.org/10.1016/j.cub.2015.11.042. 1237 24. Lu, Y.-T., Loue-Manifel, J., Bollier, N., Gadient, P., De Winter, F., Carella, P., Hoguin, 1238 A., Grey--Switzman, S., Marnas, H., Simon, F., et al. (2024). Convergent evolution 1239 of water-conducting cells in Marchantia recruited the ZHOUPI gene promoting cell 1240 wall reinforcement and programmed cell death. Curr Biol _34_, 793-807.e797. 1241 https://doi.org/10.1016/j.cub.2024.01.014. 41 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 1242 25. Lau, O.S., Song, Z., Zhou, Z., Davies, K.A., Chang, J., Yang, X., Wang, S., 1243 Lucyshyn, D., Tay, I.H.Z., Wigge, P.A., and Bergmann, D.C. (2018). Direct Control 1244 of SPEECHLESS by PIF4 in the High-Temperature Response of Stomatal 1245 Development. Curr Biol _28_, 1273-1280.e1273. 10.1016/j.cub.2018.02.054. 1246 26. Wang, L., Wan, M.C., Liao, R.Y., Xu, J., Xu, Z.G., Xue, H.C., Mai, Y.X., and Wang, 1247 J.W. (2023). The maturation and aging trajectory of Marchantia polymorpha at 1248 single-cell resolution. Dev Cell _58_, 1429-1444.e1426. 1249 10.1016/j.devcel.2023.05.014. 1250 27. Moriya, K.C., Shirakawa, M., Loue-Manifel, J., Matsuda, Y., Lu, Y.-T., Tamura, K., 1251 Oka, Y., Matsushita, T., Hara-Nishimura, I., Ingram, G., et al. (2023). Stomatal 1252 regulators are co-opted for seta development in the astomatous liverwort 1253 Marchantia polymorpha. Nature Plants _9_, 302-314. 10.1038/s41477-022-01325-5. 1254 28. Lu, K.J., van 't Wout Hofland, N., Mor, E., Mutte, S., Abrahams, P., Kato, H., 1255 Vandepoele, K., Weijers, D., and De Rybel, B. (2020). Evolution of vascular plants 1256 through redeployment of ancient developmental regulators. Proc Natl Acad Sci U 1257 S A _117_, 733-740. 10.1073/pnas.1912470117. 1258 29. Nazari, M., Fathinejad, A., Mohammadi, M.H., Bickel, S., Berg, G., and Carminati, 1259 A. (2026). The origin and evolution of root mucilage. New Phytol _249_, 1687-1698. 1260 10.1111/nph.70770. 1261 30. Ohama, N., Sato, H., Shinozaki, K., and Yamaguchi-Shinozaki, K. (2017). 1262 Transcriptional Regulatory Network of Plant Heat Stress Response. Trends Plant 1263 Sci _22_, 53-65. 10.1016/j.tplants.2016.08.015. 1264 31. Takahashi, G., Kiyosue, T., and Hirakawa, Y. (2023). Control of stem cell behavior 1265 by CLE-JINGASA signaling in the shoot apical meristem in Marchantia polymorpha. 1266 Curr Biol _33_, 5121-5131.e5126. 10.1016/j.cub.2023.10.054. 1267 32. Baumgart, L.A., Greenblum, S.I., Morales-Cruz, A., Wang, P., Zhang, Y., Yang, L., 1268 Chen, C., Dilworth, D.J., Garretson, A.C., Grosjean, N., et al. (2025). Recruitment, 1269 rewiring and deep conservation in flowering plant gene regulation. Nature Plants. 1270 10.1038/s41477-025-02047-0. 1271 33. O'Malley, R.C., Huang, S.C., Song, L., Lewsey, M.G., Bartlett, A., Nery, J.R., Galli, 1272 M., Gallavotti, A., and Ecker, J.R. (2016). Cistrome and Epicistrome Features 1273 Shape the Regulatory DNA Landscape. Cell _165_, 1280-1292. 1274 10.1016/j.cell.2016.04.038. 1275 34. Dong, Y., Krishnamoorthi, S., Tan, G.Z.H., Poh, Z.Y., and Urano, D. (2024). Co1276 option of plant gene regulatory network in nutrient responses during 1277 terrestrialization. Nature Plants _10_, 1955-1968. 10.1038/s41477-024-01851-4. 1278 35. Krishnamoorthi, S., Tan, G.Z.H., Dong, Y., Leong, R., Wu, T.-Y., and Urano, D. 1279 (2024). Hyperspectral imaging of liverwort Marchantia polymorpha identifies 1280 MpWRKY10 as a key regulator defining Foliar pigmentation patterns. Cell Rep _43_, 1281 114463. https://doi.org/10.1016/j.celrep.2024.114463. 1282 36. Cerny, M., Jedelsky, P.L., Novak, J., Schlosser, A., and Brzobohaty, B. (2014). 1283 Cytokinin modulates proteomic, transcriptomic and growth responses to 1284 temperature shocks in Arabidopsis. Plant Cell Environ _37_, 1641-1655. 1285 10.1111/pce.12270. 42 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 1286 37. Cortleven, A., Leuendorf, J.E., Frank, M., Pezzetta, D., Bolt, S., and Schmulling, T. 1287 (2019). Cytokinin action in response to abiotic and biotic stresses in plants. Plant 1288 Cell Environ _42_, 998-1018. 10.1111/pce.13494. 1289 38. Skalak, J., Cerny, M., Jedelsky, P., Dobra, J., Ge, E., Novak, J., Hronkova, M., 1290 Dobrev, P., Vankova, R., and Brzobohaty, B. (2016). Stimulation of ipt 1291 overexpression as a tool to elucidate the role of cytokinins in high temperature 1292 responses of Arabidopsis thaliana. J Exp Bot _67_, 2861-2873. 10.1093/jxb/erw129. 1293 39. Liu, B., Liu, J., Yu, J., Wang, Z., Sun, Y., Li, S., Lin, Y.J., Chiang, V.L., Li, W., and 1294 Wang, J.P. (2021). Transcriptional reprogramming of xylem cell wall biosynthesis 1295 in tension wood. Plant Physiol _186_, 250-269. 10.1093/plphys/kiab038. 1296 40. Zenker, S., Wulf, D., Meierhenrich, A., Viehöver, P., Becker, S., Eisenhut, M., 1297 Stracke, R., Weisshaar, B., and Bräutigam, A. (2025). Many transcription factor 1298 families have evolutionarily conserved binding motifs in plants. Plant Physiology 1299 _198_ . 10.1093/plphys/kiaf205. 1300 41. Yen, M.-J., Lin, K.-H., Thalimaraw, L., Yang, H.-R., Boonyaves, K., Cheng, C.-Y., 1301 and Wu, T.-Y. (2025). Conserved HSFA1-dependent chromatin dynamics drive 1302 heat stress responses in plants. Cell Rep _44_, 116714. 1303 https://doi.org/10.1016/j.celrep.2025.116714. 1304 42. Li, M., Yao, T., Galli, M., Lin, W., Zhou, Y., Chen, J.G., Gallavotti, A., and Huang, 1305 S.C. (2024). Diversification and conservation of DNA binding specificities of SPL 1306 family of transcription factors. bioRxiv. 10.1101/2024.09.13.612952. 1307 43. Strader, L., Weijers, D., and Wagner, D. (2022). Plant transcription factors — being 1308 in the right place with the right company. Current Opinion in Plant Biology _65_, 1309 102136. https://doi.org/10.1016/j.pbi.2021.102136. 1310 44. Li, M., Yao, T., Lin, W., Hinckley, W.E., Galli, M., Muchero, W., Gallavotti, A., Chen, 1311 J.-G., and Huang, S.-s.C. (2023). Double DAP-seq uncovered synergistic DNA 1312 binding of interacting bZIP transcription factors. Nat Commun _14_, 2600. 1313 10.1038/s41467-023-38096-2. 1314 45. Gremski, K., Ditta, G., and Yanofsky, M.F. (2007). The HECATE genes regulate 1315 female reproductive tract development in Arabidopsis thaliana. Development _134_, 1316 3593-3601. 10.1242/dev.011510. 1317 46. Lee, Y.S., Shiu, S.H., and Grotewold, E. (2023). Evolution and diversification of the 1318 ACT-like domain associated with plant basic helix-loop-helix transcription factors. 1319 Proc Natl Acad Sci U S A _120_, e2219469120. 10.1073/pnas.2219469120. 1320 47. Petroll, R., Varshney, D., Hiltemann, S., Finke, H., Schreiber, M., de Vries, J., and 1321 Rensing, S.A. (2025). Enhanced sensitivity of TAPscan v4 enables comprehensive 1322 analysis of streptophyte transcription factor evolution. Plant J _121_, e17184. 1323 10.1111/tpj.17184. 1324 48. Grigoriev, I.V., Hayes, R.D., Calhoun, S., Kamel, B., Wang, A., Ahrendt, S., 1325 Dusheyko, S., Nikitin, R., Mondo, Stephen J., Salamov, A., et al. (2021). 1326 PhycoCosm, a comparative algal genomics resource. Nucleic Acids Res _49_, 1327 D1004-D1011. 10.1093/nar/gkaa898. 1328 49. Feng, X., Zheng, J., Irisarri, I., Yu, H., Zheng, B., Ali, Z., de Vries, S., Keller, J., 1329 Fürst-Jansen, J.M.R., Dadras, A., et al. (2024). Genomes of multicellular algal 1330 sisters to land plants illuminate signaling network evolution. Nature Genetics _56_, 1331 1018-1031. 10.1038/s41588-024-01737-3. 43 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 1332 50. Chen, J., Luo, M., Hands, P., Rolland, V., Zhang, J., Li, Z., Outram, M., Dodds, P., 1333 and Ayliffe, M. (2023). A split GAL4 RUBY assay for visual in planta detection of 1334 protein-protein interactions. Plant J _114_, 1209-1226. 10.1111/tpj.16234. 1335 51. Bolger, A.M., Lohse, M., and Usadel, B. (2014). Trimmomatic: a flexible trimmer 1336 for Illumina sequence data. Bioinformatics _30_, 2114-2120. 1337 10.1093/bioinformatics/btu170. 1338 52. Langmead, B., and Salzberg, S.L. (2012). Fast gapped-read alignment with Bowtie 1339 2. Nat Methods _9_, 357-359. 10.1038/nmeth.1923. 1340 53. Zhang, Y., Liu, T., Meyer, C.A., Eeckhoute, J., Johnson, D.S., Bernstein, B.E., 1341 Nusbaum, C., Myers, R.M., Brown, M., Li, W., and Liu, X.S. (2008). Model-based 1342 analysis of ChIP-Seq (MACS). Genome Biol _9_, R137. 10.1186/gb-2008-9-9-r137. 1343 54. Quinlan, A.R., and Hall, I.M. (2010). BEDTools: a flexible suite of utilities for 1344 comparing genomic features. Bioinformatics _26_, 841-842. 1345 10.1093/bioinformatics/btq033. 1346 55. Bailey, T.L., Johnson, J., Grant, C.E., and Noble, W.S. (2015). The MEME Suite. 1347 Nucleic Acids Res _43_, W39-49. 10.1093/nar/gkv416. 1348 56. Hao, Y., Hao, S., Andersen-Nissen, E., Mauck, W.M., III, Zheng, S., Butler, A., Lee, 1349 M.J., Wilk, A.J., Darby, C., Zager, M., et al. (2021). Integrated analysis of 1350 multimodal single-cell data. Cell _184_, 3573-3587.e3529. 1351 10.1016/j.cell.2021.04.048. 1352 57. Hao, Y., Stuart, T., Kowalski, M.H., Choudhary, S., Hoffman, P., Hartman, A., 1353 Srivastava, A., Molla, G., Madad, S., Fernandez-Granda, C., and Satija, R. (2024). 1354 Dictionary learning for integrative, multimodal and scalable single-cell analysis. 1355 Nat Biotechnol _42_, 293-304. 10.1038/s41587-023-01767-y. 1356 58. McGinnis, C.S., Murrow, L.M., and Gartner, Z.J. (2019). DoubletFinder: Doublet 1357 Detection in Single-Cell RNA Sequencing Data Using Artificial Nearest Neighbors. 1358 Cell Syst _8_, 329-337.e324. 10.1016/j.cels.2019.03.003. 1359 59. Korsunsky, I., Millard, N., Fan, J., Slowikowski, K., Zhang, F., Wei, K., Baglaenko, 1360 Y., Brenner, M., Loh, P.-r., and Raychaudhuri, S. (2019). Fast, sensitive and 1361 accurate integration of single-cell data with Harmony. Nature Methods _16_, 12891362 1296. 10.1038/s41592-019-0619-0. 1363 60. Aran, D., Looney, A.P., Liu, L., Wu, E., Fong, V., Hsu, A., Chak, S., Naikawadi, R.P., 1364 Wolters, P.J., Abate, A.R., et al. (2019). Reference-based analysis of lung single1365 cell sequencing reveals a transitional profibrotic macrophage. Nature Immunology 1366 _20_, 163-172. 10.1038/s41590-018-0276-y. 1367 61. Wu, T., Hu, E., Xu, S., Chen, M., Guo, P., Dai, Z., Feng, T., Zhou, L., Tang, W., 1368 Zhan, L., et al. (2021). clusterProfiler 4.0: A universal enrichment tool for 1369 interpreting omics data. Innovation (Camb) _2_, 100141. 1370 10.1016/j.xinn.2021.100141. 1371 62. Phipson, B., Sim, C.B., Porrello, E.R., Hewitt, A.W., Powell, J., and Oshlack, A. 1372 (2022). propeller: testing for differences in cell type proportions in single cell data. 1373 Bioinformatics _38_, 4720-4726. 10.1093/bioinformatics/btac582. 1374 63. Cao, J., Spielmann, M., Qiu, X., Huang, X., Ibrahim, D.M., Hill, A.J., Zhang, F., 1375 Mundlos, S., Christiansen, L., Steemers, F.J., et al. (2019). The single-cell 1376 transcriptional landscape of mammalian organogenesis. Nature _566_, 496-502. 1377 10.1038/s41586-019-0969-x. 44 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 1378 64. Qiu, X., Mao, Q., Tang, Y., Wang, L., Chawla, R., Pliner, H.A., and Trapnell, C. 1379 (2017). Reversed graph embedding resolves complex single-cell trajectories. 1380 Nature Methods _14_, 979-982. 10.1038/nmeth.4402. 1381 65. Trapnell, C., Cacchiarelli, D., Grimsby, J., Pokharel, P., Li, S., Morse, M., Lennon, 1382 N.J., Livak, K.J., Mikkelsen, T.S., and Rinn, J.L. (2014). The dynamics and 1383 regulators of cell fate decisions are revealed by pseudotemporal ordering of single 1384 cells. Nat Biotechnol _32_, 381-386. 10.1038/nbt.2859. 1385 66. Street, K., Risso, D., Fletcher, R.B., Das, D., Ngai, J., Yosef, N., Purdom, E., and 1386 Dudoit, S. (2018). Slingshot: cell lineage and pseudotime inference for single-cell 1387 transcriptomics. BMC Genomics _19_, 477. 10.1186/s12864-018-4772-0. 1388 67. Van den Berge, K., Roux de Bézieux, H., Street, K., Saelens, W., Cannoodt, R., 1389 Saeys, Y., Dudoit, S., and Clement, L. (2020). Trajectory-based differential 1390 expression analysis for single-cell sequencing data. Nat Commun _11_, 1201. 1391 10.1038/s41467-020-14766-3. 1392 68. Aibar, S., González-Blas, C.B., Moerman, T., Huynh-Thu, V.A., Imrichova, H., 1393 Hulselmans, G., Rambow, F., Marine, J.-C., Geurts, P., Aerts, J., et al. (2017). 1394 SCENIC: single-cell regulatory network inference and clustering. Nature Methods 1395 _14_, 1083-1086. 10.1038/nmeth.4463. 1396 69. Gu, Z., Eils, R., and Schlesner, M. (2016). Complex heatmaps reveal patterns and 1397 correlations in multidimensional genomic data. Bioinformatics _32_, 2847-2849. 1398 10.1093/bioinformatics/btw313. 1399 70. Chen, S., Zhou, Y., Chen, Y., and Gu, J. (2018). fastp: an ultra-fast all-in-one 1400 FASTQ preprocessor. Bioinformatics _34_, i884-i890. 10.1093/bioinformatics/bty560. 1401 71. Love, M.I., Huber, W., and Anders, S. (2014). Moderated estimation of fold change 1402 and dispersion for RNA-seq data with DESeq2. Genome Biol _15_, 550. 1403 10.1186/s13059-014-0550-8. 1404 72. Michna, A., Braselmann, H., Selmansberger, M., Dietz, A., Hess, J., Gomolka, M., 1405 Hornhardt, S., Blüthgen, N., Zitzelsberger, H., and Unger, K. (2016). Natural Cubic 1406 Spline Regression Modeling Followed by Dynamic Network Reconstruction for the 1407 Identification of Radiation-Sensitivity Gene Association Networks from Time1408 Course Transcriptome Data. PLoS One _11_, e0160791. 1409 10.1371/journal.pone.0160791. 1410 73. Supek, F., Bošnjak, M., Škunca, N., and Šmuc, T. (2011). REVIGO summarizes 1411 and visualizes long lists of gene ontology terms. PLoS One _6_, e21800. 1412 10.1371/journal.pone.0021800. 1413 74. Morabito, S., Reese, F., Rahimzadeh, N., Miyoshi, E., and Swarup, V. (2023). 1414 hdWGCNA identifies co-expression networks in high-dimensional transcriptomics 1415 data. Cell Reports Methods _3_, 100498. 1416 https://doi.org/10.1016/j.crmeth.2023.100498. 1417 75. Shahan, R., Hsu, C.W., Nolan, T.M., Cole, B.J., Taylor, I.W., Greenstreet, L., Zhang, 1418 S., Afanassiev, A., Vlot, A.H.C., Schiebinger, G., et al. (2022). A single-cell 1419 Arabidopsis root atlas reveals developmental trajectories in wild-type and cell 1420 identity mutants. Dev Cell _57_, 543-560 e549. 10.1016/j.devcel.2022.01.008. 1421 76. Chang, G., Ma, J., Wang, S., Tang, M., Zhang, B., Ma, Y., Li, L., Sun, G., Dong, S., 1422 Liu, Y., et al. (2023). Liverwort bHLH transcription factors and the origin of stomata 45 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 1423 in plants. Curr Biol _33_, 2806-2813.e2806. 1424 https://doi.org/10.1016/j.cub.2023.05.050. 1425 1426 46 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 1427 **Figures** 47 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 1428 48 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 1429 49 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 1430 50 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 1431 51 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 1432 52 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 1433 53 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 1434 **Supporting Figures** 1435 54 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 1436 55 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 1437 56 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 1438 57 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 1439 58 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 1440 59 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 1441 60 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 1442 61 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 1443 62 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 1444 63 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 1445 64 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 1446 65 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 1447 66 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 1448 67 [bioRxiv preprint doi: https://doi.org/10.64898/2026.11.09.765476; this version posted November 12, 2026. The copyright holder for this preprint](https://doi.org/10.64898/2026.11.09.765476) (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made [available under a CC-BY-NC 4.0 International license.](http://creativecommons.org/licenses/by-nc/4.0/) 1449 68