Diversification of heat shock transcription factors expanded thermal stress responses during early plant evolution.

The copy numbers of many plant transcription factor (TF) genes substantially increased during terrestrialization. This allowed TFs to acquire new specificities and thus create gene regulatory networks (GRNs) with new biological functions to help plants adapt to terrestrial environments. Through characterizing heat shock factor (HSF) genes MpHSFA1 and MpHSFB1 in the liverwort Marchantia polymorpha, we explored how heat-responsive GRNs widened their functions in M. polymorpha and Arabidopsis thaliana. An interspecies comparison of heat-induced transcriptomes and the evolutionary rates of HSFs demonstrated the emergence and subsequent rapid evolution of HSFB prior to terrestrialization. Transcriptome and metabolome analyses of M. polymorpha HSF-null mutants revealed that MpHSFA1 controls canonical heat responses such as thermotolerance and metabolic changes. MpHSFB1 also plays essential roles in heat responses, as well as regulating developmental processes including meristem branching and antheridiophore formation. Analysis of cis-regulatory elements revealed development- and stress-related TFs that function directly or indirectly downstream of HSFB. Male gametophytes of M. polymorpha showed higher levels of thermotolerance than female gametophytes, which could be explained by different expression levels of MpHSFA1U and MpHSFA1V on sex chromosome. We propose that the diversification of HSFs is linked to the expansion of HS responses, which enabled coordinated multicellular reactions in land plants.

Introduction

Simple gene regulatory mechanisms in ancestral organisms have expanded to higher complexity as species evolved. Such advancement in gene regulatory systems often occurs in a saltatory (discontinuous) manner, allowing for the rapid acquisition of new biological roles. Transcription factors (TFs) are crucial for the transcriptional regulation of numerous genes. Therefore, the numbers and combinations of TF genes expressed in cells contribute to the diversity of gene regulatory networks (GRNs) (Chalancon et al., 2012; Marshall-Colón and Kliebenstein, 2019). Cis-regulatory elements determine the architectures and functions of GRNs, which gradually evolved during early plant life (Erwin and Davidson, 2009; Stergachis et al., 2014; Wu et al., 2021).

Heat shock transcription factors (HSFs) are central regulators of heat responses such as protein biosynthesis and folding and metabolic changes in all eukaryotes. In addition to these conventional functions, plant HSFs have acquired broader roles in growth regulation and biotic stress responses (Scharf et al., 2012), possibly through gene duplications and subsequent mutations in angiosperms (Lehti-Shiu et al., 2017). Reflecting their broader roles, land plants contain multiple HSF genes (e.g. 21 HSF genes in Arabidopsis thaliana) that exhibit functional distinctions between knockout phenotypes (Scharf et al., 2012). These distinct functions can be explained by tissue- and stress type-specific induction mechanisms as well as their diversified domain structures, which classify HSFs into three families: HSFA, B, and C (Andrási et al., 2021). Members of the HSFA subclass serve as master transcriptional activators of heat stress (HS) responses. Genetic ablation of Arabidopsis HSFA1 caused impaired heat-sensitive phenotypes and altered gene expression (Yoshida et al., 2011; Liu et al., 2011). In contrast, HSFB family member act as either co-activators repressors of HSFA, depending on the species and/or experimental conditions (Bharti et al., 2004; Ikeda et al., 2011a; Fragkostefanakis et al., 2019a). In addition to HS responses, the ectopic expression of HSFB genes alters root development under normal growth conditions (Begum et al., 2013).

Morphological control accounts for the developmental adaptation strategies of land plants under HS or high ambient temperature. These unfavorable thermal conditions lead to elongated stems and hyponasty in some crops, in turn leading to a reduction in total biomass (Hasanuzzaman et al., 2013). HS also induces the accumulation of certain metabolites as protectants against high temperatures. More specifically, HS alters the core metabolic fluxes of glycolysis, the pentose phosphate pathway, the TCA cycle, and starch and amino acid biosynthesis pathways (Raza, 2020). Nevertheless, a thorough understanding of metabolic gene expression and metabolite accumulation under HS conditions and the involvement of HSFs is yet to be achieved, mainly due to the high gene redundancy in flowering plants (Lehti-Shiu et al., 2017).

Here, using interspecies meta-transcriptome analysis and HSF knockout mutants in Marchantia polymorpha, we analyzed the diversification and expansion of heat-responsive GRNs during aquatic to land plant evolution. Phylogenetic and evolutionary rate analyses of HSFs revealed the emergence of HSF subclasses at the base of Streptophyta and the strong accumulation of mutations in HSFB genes. Marchantia polymorpha HSFA1 and HSFB1 knockout mutants (Mphsfa1 and Mphsfb1) showed distinct transcriptomic defects under both HS and non-HS conditions, which can be explained by the different usage of cis-regulatory sequences from direct and/or indirect target genes. The Mphsfa1 mutants showed increased thermal sensitivity, such as thallus senescence, metabolic alterations, and growth retardation. The Mphsfb1 mutants showed altered expression profiles of genes related to developmental processes as well as typical heat responses. The unconventional role of MpHSFB1 was confirmed by examining a knockout mutant that displayed developmental defects in meristem and reproductive organ formation. We also uncovered gender- and developmental stage-based differences in HS response in M. polymorpha, which could be explained by female- and male-specific MpHSFA1 genes on U and V chromosomes (MpHSFA1U and MpHSFA1V, respectively). Collectively, the diversification of HSFs in early land plants has expanded their gene regulatory and biological roles associated with developmental and metabolic responses to HS.

Results

The functions of HS-responsive genes expanded during early plant evolution

It is widely accepted that land plants acquired greatly advanced morphological and intracellular mechanisms to adapt to temperature fluctuations. To compare the HS-responsive mechanisms from green algae to angiosperms, we first analyzed public and in-house transcriptome datasets from a Streptophyta alga (Klebsormidium nitens), a bryophytic plant (M. polymorpha), three flowering plants (A. thaliana, tomato [Solanum lycopersicum] and rice [Oryza sativa]), and the yeast Saccharomyces cerevisiae (Supplemental Figure S1A and Supplemental Data Set 1). These species represent different taxa distributed along the plant evolutionary timescale. K. nitens is a close relative to land plants (Lewis and McCourt, 2004; de Vries and Archibald, 2018), M. polymorpha represents a bryophytic plant, and A. thaliana, S. lycopersicum, and O. sativa are dicot or monocot models that have been extensively studied (Kotak et al., 2007; Mittler et al., 2012).

The multi-taxa analysis yielded approximately 700 to approximately 2,600 HS-responsive genes in different species (Supplemental Figure S1A and Supplemental Data Set 1). Orthologous gene analysis of the HS-responsive genes (Supplemental Figure S1B and S1C) identified 10 and 70 HS-responsive orthogroups conserved in all six species and in three angiosperms, respectively (Figure 1, A and B, highlighted in red text). A higher overlap of HS-induced gene orthologs was observed within four land plants (the degree of overlap ranged from 0.14 to 0.3), which was greater than the median of the degree of overlap between all pairs of six species (median = 0.125), while HS-repressed genes tended to be conserved from K. nitens to angiosperms (Figure 1, C and D). HS-induced genes conserved in all six species and in three plants were associated with the gene ontology (GO) terms abiotic stress response and metabolic process (magenta text, Figure 1E) and regulation of gene expression, development, or reproduction (green text, Figure 1E), respectively. On the other hand, HS-repressed genes tended to be associated with metabolic process and cell growth in all species examined (Supplemental Figure S1D).

Figure 1

Transcriptomic comparison reveals shock (HS)-responsive genes that were conserved in Eukaryotes and acquired in green plants. A and B, UpSet plots showing the number of HS-induced (A) or HS-repressed (B) orthogroups in S. cerevisiae (yeast), K. nitens (green alga), M. polymorpha (liverwort), A. thaliana (Arabidopsis), S. lycopersicum (tomato), and O. sativa (rice). C and D, Pairwise comparison of the degree of overlap of HS-induced (C) or HS-repressed (D) orthogroups among the selected taxa. Colors in the heatmaps represent the degree of overlap (red: higher; white: lower). Green boxes indicate that the degree of overlap is greater than the median value (C: 0.125; D: 0.157). E, GO terms (Biological process, BP) identified using Arabidopsis orthologs in HS-induced overlapping orthogroups between Arabidopsis and other species, with red indicating greater enrichment (−log10(FDR)). The numbers in each cell represent the numbers of Arabidopsis orthologs associated with GO terms shown on the right. Magenta text indicates GO terms identified in all species analyzed, while green text indicates GO terms identified in land plant species.

Given that HSFs evolved to form sub-families in land plants (Scharf et al., 2012; Wang et al., 2018) and that HSFs are the master regulators of HS responses in plants (Ohama et al., 2017), we hypothesized that the divergence of HS-induced genes was associated with the fixation and specialization of HSFs in plants. We observed that HSFA1 has maintained the original residues more strictly (turquoise, Figure 2A), while HSFB1 spans a large phylogenetic distance from its emergence to the separation of these proteins in Streptophyta species (Figure 2, A and B; P < 0.01). Analysis of the deduced sequences of ancestral HSFs (Supplemental Figure S2A) revealed that 3 and 10 residues in the HSFA and HSFB subclasses, respectively, mutated from the common HSF ancestor. The mutated residues were located in the α-helix 2 (H2) in the DNA binding domain (DBD), where the domain interacts with DNA (Supplemental Figure S2B). In general, members of the HSFA and HSFB subclasses contain conserved surface residues for DNA-binding in land plants (Supplemental Figure S2C), although the H2 helix in HSFA underwent small but further evolutionary changes (Supplemental Figure S2D). The HSFC subclass was separated from HSFA in flowering plants, showing a large branch distance within the HSFA subtree (light green, Figure 2A). These results suggest that the plant HSFB and HSFC subclasses accumulated multiple mutations during the transition from aquatic to land plants and within angiosperms, which might be related to the increasing complexity of plant HS responses.

Figure 2

HSFA1 and HSFB1 have different phylogenetic distances in plants. A, Phylogeny of plant HSF sequences. The maximum-likelihood tree was constructed using the DBD sequences of selected taxa. HSFA, HSFB, and HSFC groups are shown in different colors. Red text denotes MpHSFA1U, MpHSFA1V and MpHSFB1. Colored nodes indicate HSF sequences in A. thaliana (red), O. sativa (green), bryophyta (purple), and streptophycean algae (yellow). All sequences used to construct the phylogenetic trees are shown in Supplemental File 1 and 2. B, Duplication length from the HSFA and HSFB groups. **P < 0.01, determined by Mann–Whitney U Test. C, Heatmap (left panel) showing the enrichment (−log10(FDR)) of HS-induced TF families in five species. Colors in the heatmap represent the significance determined by hypergeometric distribution followed by BH adjustment. Heat map showing the enrichment (−log10(FDR)) of predictive TF binding motifs determined by Wilcoxon signed rank test using AME suites. Gray cells denote not applicable (NA); red text represents significantly enriched TFs and motifs found in at least four species.

The number of TFs gradually increased from aquatic to terrestrial plants (Bowman et al., 2017; Wilhelmsson et al., 2017) (Supplemental Figure S2E). Gene number comparisons in aquatic and land plants revealed 40 TF families that are present in at least one green alga genome (Supplemental Figure S2E). Of these, 22 TF families including HSF multiplied from one family in green algae into two or more sub-classes in land plants (Supplemental Figure S2E). To evaluate the conservation and diversification of their roles in HS, we performed enrichment analysis of TF families in the meta-transcriptome datasets. Although several TF families exhibited a higher chance of being induced/repressed under HS in at least one species (Figure 2C, left), only the HSF family showed a conserved HS response in all five species except for S. lycopersicum (Figure 2C, left). The HS-regulated TFs tended to contain cis-regulatory elements for HSFs and other TFs (Figure 2C, right). These results highlight the central role of HSFs in regulating other TFs either directly or indirectly under HS and the HS-dependent regulation of stress-responsive GRNs.

MpHSFs regulate distinct developmental processes and thermotolerance

Marchantia polymorpha are bryophytic plants that possess HSFA and HSFB (Wang et al., 2018) (Figure 2A, red text), and MpHSFA1 and MpHSFB1 are ubiquitously expressed in most tissues in Marchantia (Supplemental Figure S2F). To elucidate the functional differences between the two HSF subfamilies, we generated CRISPR/Cas9 knockout mutants for the single copy MpHSF genes (Mphsfa1V and Mphsfb1 mutants) with different indel editing sites (Supplemental Figure S2G). We first analyzed the transcriptomic and metabolic profiles of the mutants under control and HS conditions (Figures 3 and 4); the morphological phenotypes and HS responses of the mutants are described in Figures 5–7.

Figure 3

HSFA1 mainly regulates canonical heat-responsive genes, and HSFB1 regulates genes involved in developmental processes and defense responses. A, Overview of HS-responsive genes regulated by MpHSFA1V and MpHSFB1. Transcriptomic profiles were analyzed in 7-day-old Mphsfa1V, Mphsfb1, and Tak-1 plants treated at 37°C for 1 h. UpSet plots show the number of up- and down-regulated genes (cut off threshold, |log2 Fold Change (log2(FC))| ≥ 1; FDR < 0.1), demonstrating the different expression patterns of Mphsfa1V-specific and Mphsfb1-psecific genes (top bar graphs). The total number of up- and down-regulated genes in each mutant is shown on the right. B, Expression patterns of common HS-responsive DEGs in Mphsfa1V, Mphsfb1, and Tak-1 plants. Black dots represent the fold increase (top violin plot) or decrease (bottom) of common HS-responsive genes between control and HS conditions. Red bars represent the median. ***P < 0.001, determined by Mann–Whitney U Test, compared to Tak-1. C, GO terms associated with the HS response and/or Mphsfa1V or Mphsfb1 null mutation. Color in the heatmap indicates the level of enrichment (−log10(FDR)). Genes whose expression changes were four-fold smaller in Mphsfa1V or Mphsfb1 than in Tak-1 plants under HS were classified into Mphsfa1V- and Mphsfb1-dependent genes (red and blue dots represent up- and downregulation, respectively). GO terms related to metabolism changes, developmental changes, and the HS response are highlighted in blue, green, and magenta, respectively. D and E, Treemaps summarizing GO terms enriched in Mphsfa1V-specific genes (D.) or Mphsfb1-specific genes (E.) from C (include up- and down-regulated genes). Box sizes in the treemaps represent the degree of significance (FDR). Red text highlights the GO terms that were enriched only in either Mphsfa1V- or Mphsfb1-specific genes. F, Enrichment of orthologous Mphsfa1V-regulated genes in M. polymorpha (this study) and A. thaliana (dataset reanalyzed from (Liu et al., 2011; Liu and Charng, 2013), Supplemental Figure S1A). P-value was determined by Fisher’s chi-square test. G, Top 10 GO BP terms of overlapping genes from (F.).

Figure 4

Heat-induced changes in primary metabolism differ between Mphsfa1 and Mphsfb1. A, Thermal regulation of genes involved in primary metabolism in wild type and Mphsf plants. Dot plots show the median log2FC of heat-dependent induction (red plots) or suppression (blue plots) of metabolite-related genes in Mphsfa1V, Mphsfb1, and Tak-1 plants. Metabolite-related genes were categorized using the PlantCyc database. PlantCyc pathways with three or more DEGs in at least one genotype are shown. *P < 0.05, **P < 0.01, and ***P < 0.001, determined by Mann–Whitney U Test compared to Tak-1 plants. Magenta text represents genes whose expression was mis-regulated in both Mphsfa1V and Mphsfb1, while green text indicates that the difference was observed only in Mphsfa1V. B, Dynamics of metabolic changes in Tak-1, Mphsfa1V, and Mphsfb1 plants under HS treatment at 37°C for 1 and 5 h. The heatmap shows the relative changes in sugar, sugar alcohol, amino acid, and organic acid contents determined by GC–MS. Colors indicate FCs in the Z-Score values of the relative metabolite contents across all tested conditions (n = 3; *P < 0.05, as determined by two-tailed Student’s t test).

Figure 5

Mphsfa1V and Mphsfb1 mutations caused distinct morphological defects and thermal tolerance. A and B, Representative photographs of Mphsfa1v (A) (Mphsfa1V-1, Mphsfa1V-13 and Mphsfa1V-19) and Mphsfb1 (B) (Mphsfb1-25, Mphsfb1-29 and Mphsfb1-70) lines compared to Tak-1 plants. Genome editing sites for the three independent lines are shown in Supplemental Figure S2G. The photographs were taken 25 days after sowing gemma on 1/2 × B5 plates at 22°C. C and D, Thallus area of Mphsfa1V (C) and Mphsfb1 (D) mutants growing on 1/2 × B5 plates over time. Dots and curves represent data from individual plant and the fitted model using local regression analysis, respectively (n = 5–10 independent plants per line). ***P < 0.001, determined by Kruskal–Wallis test and Mann–Whitney U test with Bonferroni correction. E and F, Thermal tolerance phenotypes of Tak-1, Mphsfa1V, and Mphsfb1 plants. Plant images were taken after HS treatment at 37°C for 24 h with a subsequent recovery period of 1, 14, and/or 21 days at 22°C. Note that plants in (E) and (F) were derived from the same plates but photographed after different recovery periods (14 versus 21 Days at 22°C). G, Ratio of green to total thallus area of three Mphsfa1V lines and Tak-1 plants on Days 1 and 14 are shown with boxplots. Boxplots present all data (black dots) with the median values, the first and third quartiles, and whiskers of maximum and minimum values. H, Thermal tolerance phenotypes of Mphsfb1 mutants following 2-h heat shock treatment. Plant images were taken after HS treatment at 37°C for 2 h in a water bath and recovery at 22°C for 5 days. I, Survival rates of the Mphsfb1 and Tak-1 plants are presented in boxplots. Boxplots present all data (black dots) with the median values, the first and third quartiles and whiskers of maximum and minimum values. The experiments were repeated at least four times, and 20 gemmalings from each line were tested each time. *P < 0.05; **P < 0.01, and ***P < 0.001, determined by the two-tailed Student’s t test. Morphology of the Mphsfa1Vb1 double knockout mutant is shown in Supplemental Figure S6, E–G.

Figure 6

Female Mphsfa1U mutants exhibited growth defects and high HS sensitivity. A, Morphological defects in Mphsfa1U (Mphsfa1U-13, Mphsfa1U-21, and Mphsfa1U-29) and BC3 plants under non-HS conditions. Genome editing sites of the mutant lines are presented in Supplemental Figure S7C. Images show representative plants grown on 1/2 × B5 plates for 25 days at 22°C. B, Thermal tolerance of BC3 and Mphsfa1U plants. Images show representative plants treated at 37°C for 24 h in a growth chamber with a subsequent recovery at 22°C for 1 or 14 days. The quantitative data are presented in Figure 6C (right boxplots). C, Thallus area of BC3 and Mphsfa1U plants growing on 1/2 × B5 plates over time (left). The graph shows the thallus area of 5–10 individual plants per genotype (dots) with a representative growth curve. **P < 0.01, determined by Kruskal–Wallis test and Mann–Whitney U test with Bonferroni correction. Right boxplots show the ratio of green area to total thallus area on Days 1 and 14 after HS treatment. Boxplots present all data (black dots) with the median values, the first and third quartiles, and whiskers of maximum and minimum values. The experiments were repeated at least four times, with 20 gemmalings per genotype in each experiment. *P < 0.05; **P < 0.01, and ***P < 0.001, determined by two-tailed Student’s t test. D, Global similarity of transcriptome from each genotype under control (Cntr) and HS conditions. DEGs identified in at least one genotype or condition were used for PCA. The first two components are shown. E and F, HSFA1-dependent heat-responsive genes in Tak-1 or BC3. Intersect of the Venn diagram indicates the number of heat-responsive genes that showed similar expression changes due to the Mphsfa1V and Mphsfa1U mutations in both genders (within two-fold differences between Tak-1/Mphsfa1V and BC3/Mphsfa1U). GO terms enriched in the commonly regulated genes (n = 312, red text) are shown in (F). G and H, Global similarity of metabolome changes in wild type and mutant lines under control (Cntr) and HS conditions. Metabolites that differentially accumulated in at least one genotype or condition were used for PCA. The first and third components are shown in G. Factors contributing to the first and third components are shown in (H). Metabolites commonly regulated in Mphsfa1V and Mphsfa1U plants (n = 25, red text in Venn diagram) were used for the factor analysis.

Transcriptome analysis revealed over 3,000 differentially expressed genes (DEGs) in Mphsfa1V and/or Mphsfb1 under HS (Figure 3A, Supplemental Figure S3A, Supplemental Data Set 2). Under control conditions (Supplemental Figure S3B), 40 genes (6.9%) overlapped among Mphsfa1V- and Mphsfb1-specific genes, while 479 genes (82.6%) were specifically mis-regulated in Mphsfb1. A large portion of HS-responsive genes detected in the wild type (Figure 1) still exhibited HS responses in Mphsfa1V and Mphsfb1 (blue columns, Figure 3A), although their induction or repression levels were smaller in the mutants (Figure 3B; P < 0.001). This suggests the importance of both MpHSF genes and the existence of an HSF-independent mechanism during thermal responses. Among the 1,124 HS-regulated DEGs in Mphsfb1 and 1,076 in Mphsfa1V, 473 (42%) and 300 (27.9%) HS-regulated DEGs overlapped with HS-responsive DEGs in wild-type Tak-1, respectively (Supplemental Figure S3B). These results suggest that although MpHSFA1 and MpHSFB1 have individual specialized functions, their functions in regulating HS responses are partially redundant.

GO enrichment analysis revealed that MpHSFA1V regulates canonical HS responses, such as the unfolded protein response, chaperone, and hormone signaling (Figure 3C, magenta text and D), whereas MpHSFB1’s functions are related to cell division, cell-wall modification, and transporters (Figure 3C, green text and E, false discovery rate (FDR) cutoff = 0.1). The requirements of HSF for both basal and HS-regulated gene expression were further confirmed by qRT-PCR of MpHSP90, MpHSP70, MpsHSP, MpHSFA1V, and MpHSFB1 (Supplemental Figure S3, B and C). The roles of HSFA1 in transcriptional regulation appear to be conserved in M. polymorpha and A. thaliana, as shown by the significant overlap of the MpHSFA1V-specific genes with those of the Arabidopsis hsfa1 mutant (P < 1 × 10−7; Figure 3, F and G).

Promoter analysis (within 2-kb upstream of the transcription starting site) of these HS-responsive genes revealed the presence of binding motifs for AP2EREBP, C2H2, C2C2-dof, and BBRBPC TFs (red text in Supplemental Figure S3D), with a higher probability than expected. Interestingly, heat shock elements (HSEs) were significantly enriched only in the promoter regions of Mphsfa1V-specific genes (33% of Mphsfa1V-specific genes) under non-HS conditions, while cis-regulatory motifs for Golden2 (G2)-like TFs and WRKY TFs were more highly enriched in Mphsfb1-specific genes (Supplemental Figure S3D). We then analyzed cis-regulatory elements in the promoter regions of Mphsfa1V-specific or Mphsfb1-specific TF genes (Supplemental Figure S3E). HSE was predominantly enriched in Mphsfa1V-specific TF genes, while NAM, ATAF and CUC (NAC), basic Helix-Loop-Helix (bHLH), and Teosinte branched1/Cincinnata/proliferating cell factor (TCP)-binding motifs were more significantly found in Mphsfb1-dependent TF genes.

Transcriptional regulation occurs for both direct and indirect targets. To further understand the functions and motif sequence patterns of genes that are likely the direct target of MpHSFA1 and MpHSFB1, we narrowed the analysis to DEGs possessing HSE in their promoter regions. First, we analyzed DEGs downregulated by Mphsfa1V and Mphsfb1 mutations under control conditions (22 and 202 DEGs, respectively; Supplemental Figure S4A). Of these DEGs, 11 DEGs in Mphsfa1V (11/22 = 50%) and 70 DEGs in Mphsfb1 (70/202 = 34.7%) contained at least one HSE in their promoter regions (Supplemental Figure S4A). GO analysis revealed that the 11 Mphsfa1V-specific DEGs were associated with the heat response, abiotic stress response, and protein folding, while the 70 Mphsfb1-specifc DEGs were related to ion transport, homeostatic process, and morphogenesis (FDR < 0.05, Supplemental Figure S4B).

We then extracted HS-induced DEGs containing at least one HSE (Supplemental Figure S4C), as HSEs are generally required for the activation of HS-induced genes (Supplemental Figure S4C). The analysis identified HSE-containing DEGs: 205 genes (205/332 = 61.7%) mis-regulated specifically in Mphsfa1V and 56 genes (56/130 = 43.1%) in Mphsfb1, where the lower percentage in Mphsfb1 suggested more indirect regulation by MpHSFB1 (Supplemental Figure S4C). GO analysis revealed that DEGs not induced by HS in Mphsfa1V were associated with canonical HS responses, such as protein folding and chaperone activities, while DEGs not induced in Mphsfb1 were associated with the wounding response, metabolic process, and auxin transport (FDR < 0.05, Supplemental Figure S4D). We next analyzed HSE motif sequences in the Mphsfa1V- or Mphsfb1-specific DEGs, which would represent HSFA1-binding and HSFB1-binding motifs, respectively (Supplemental Figure S4, E and F). Twelve nonredundant HSEs were identified as HSFA1-binding motifs, nine of which contained a canonical HSE sequence (GAAnnTTC) (Supplemental Figure S4E). On the other hand, six nonredundant HSEs were identified as HSFB1-binding motifs, but only two sequences possessed the canonical HSE sequence (Supplemental Figure S4E). These results suggest that MpHSFA1 regulates canonical HS responses mainly through direct target genes, while MpHSFB1 controls specialized HS responses both from direct binding (wounding, metabolic process, and auxin transport) and indirect binding (cell wall biosynthesis and development) under HS conditions.

Analysis of the top 10 enriched motifs from MpHSFA1V overexpression chromatin immunoprecipitation (ChIP)-seq peaks validated HSF-binding elements as the most significantly enriched motifs (Supplemental Figure S5, A and B, FDR = 1.01 × 10−11). A few other motifs similar to ERF, bZIP, C2C2, and BBRBPC binding motifs were also identified (FDR = 5.08 × 10−11 ∼ 1.40 × 10−7, Supplemental Figure S5A). Furthermore, mis-regulated genes in the Mphsfa1V mutants significantly overlapped with genes directly bound by HSFA1 (117, FDR = 1.15 × 10−11), suggesting they are direct regulatory targets of HSFA1 (Supplemental Figure S5C). GO analysis of the direct regulatory targets revealed functions associated with protein folding, response to stress, and response to HS (FDR < 0.05 as cut off, Supplemental Figure S5D). Collectively, these results possibly attribute both functional overlaps and distinctions between MpHSFA1V and MpHSFB1 to the binding preference of cis-regulatory motifs from direct or indirect target genes and TFs. These motifs are associated with the overlap of either the HSFA or HSFB subclass (Supplemental Figures S3E and S4) with the motifs of heat-induced TFs across five plant species (Figure 2C). This suggests that each HSF subclass acquired specialized functions by recruiting a different set of TFs in their GRNs to trigger HS-responsive cascades either commonly conserved in multiple species or only in specific plant lineages.

MpHSFs modulate primary metabolite production in response to HS

HS modulates metabolic processes in various eukaryotes (Figures 1 and 2), which are likely mediated by HSFs in land plants (Figure 3). To further dissect the conserved HS responses, we extracted metabolic pathway genes from the PlantCyc database (Schläpfer et al., 2017; Figure 4A). Compared to Tak-1 plants, Mphsfa1V and Mphsfb1 plants exhibited a reduced number and extent of HS-induced or repressed metabolic genes involved in GO categories amino acid, hormone, secondary metabolite, fatty acid, nucleoside and nucleotide, and transport protein biosynthesis or degradation (Figure 4A, magenta text; P < 0.05). Extraction of PlantCyc metabolic genes from meta-transcriptome analysis revealed that similar metabolic adjustments under HS are conserved from bryophytic to flowering plants (Supplemental Figure S5E). Interestingly, the expression of metabolic genes involved in energy pathways such as photosynthesis and aerobic respiration were drastically altered in Mphsfa1 but not in Mphsfb1 plants (Figure 4A, green text), highlighting additional functional distinctions between the HSFA and HSFB subclasses.

These transcriptome-based inferences were confirmed by profiling primary metabolites by gas chromatography–mass spectrometry (GC–MS). The analysis yielded nearly 300 peaks regulated by HS, which contained 80 compounds annotated in NIST mass spectral library (Figure 4B;Supplemental Data Set 3). In line with the transcriptomic changes in anabolic pathway genes, HS resulted in higher accumulation of various sugars and amino acids such as alanine, proline, phenylalanine, and tyrosine (Figure 4B; P < 0.05, compared to non-HS control conditions). Compared to Tak-1, Mphsfa1V and Mphsfb1 plants exhibited higher steady-state levels of these amino acids under non-HS conditions, with lower levels of induction under HS (Figure 4B, ∼70% and ∼39% of measured amino acids that significantly accumulated in Tak-1 versus Mphsfa1V or Mphsfb1 plants, respectively). In contrast, in Mphsfa1V and Mphsfb1 plants, ∼40% of the measured organic acids highly accumulated after 5 h of HS treatment. Of note, HS induced GABA accumulation only in Tak-1 and Mphsfb1 but not in Mphsfa1V, suggesting that MpHSFA1 plays unique roles in the metabolic HS response (Figure 4B; P < 0.05).

Metabolic adjustments are often linked to developmental strategies under stress conditions. Therefore, we quantified morphological changes and growth retardation in the Mphsf mutants (Figure 5, A–D). Under non-HS conditions, Mphsfa1V (Figure 5A) and Mphsfb1 (Figure 5B) plants exhibited a 1.5-fold (Figure 5C) and 4-fold (Figure 5D) smaller thallus area than the wild-type Tak-1 (P < 0.001), with early and slightly delayed senescence of the thallus, respectively (Supplemental Figure S6, A and C). In addition, the Mphsfb1 mutants, but not Mphsfa1V, showed a higher rhizoid-to-thallus ratio, a general indicator of stress in M. polymorpha (Supplemental Figure S6, B and D; P < 0.001). A thermotolerance assay in a 37°C growth chamber for 24 h revealed increased HS sensitivity of Mphsf1V (Figure 5, E–G; P < 0.001) but not Mphsfb1 (Figure 5F). Nevertheless, Mphsfb1 mutants exhibited hypersensitivity to HS when they were subjected to an acute 37°C treatment for 2 h in a water bath (Figure 5, H and I). These results suggest that MpHSFA1 and MpHSFB1 are both required for thermotolerance but have distinct roles in morphological stress responses.

Two MpHSFA1 genes encoded on the U and V chromosomes have comparable functions in thallus development and HS responses

MpHSFA1 is encoded as a single-copy gene on the U and V chromosomes of female and male gametophytes (BC3 and Tak-1, respectively) (Bowman et al., 2017). We hypothesized that there are gender-based differences in thermal responses associated with the female- or male-specific HSFA1. A whole transcriptome comparison between BC3 and Tak-1 revealed more than 1,000 HS-responsive genes shared by the two gametophytes (Supplemental Figure S7, A and B and Supplemental Data Set 4). The HS-induced genes that overlapped between BC3 and Tak-1 were related to the GO terms heat response, light response, abiotic stress, and protein folding (All, Supplemental Figure S7B). In both gametophytes, temperature response genes showed acute expression changes within 1 h of HS, whereas the expression of metabolism process genes was altered later, at 5 h of HS treatment (Supplemental Figure S7B, FDR cut-off: < 0.1). These results indicate that female and male gametophytes have comparable transcriptomic and metabolic responses to HS.

The functions of MpHSFA1U were genetically confirmed by generating female Mphsfa1U mutants (Figure 6 and Supplemental Figure S7, C–F). Compared to the parental line BC3, the Mphsfa1U mutants exhibited growth retardation and a smaller green area (Figure 6, A and C, left; P < 0.01) and early senescence (Supplemental Figure S7, D and F) but no change in rhizoid coverage (Supplemental Figure S7, E and G). Resembling Mphsfa1V mutants (Figure 3), the Mphsfa1U mutants showed reduced HS tolerance. After 37°C treatment for 24-h with a recovery period of 14 days, 75%–100% of the thallus area remained green in BC3 female plants compared to 0%–50% in Mphsfa1U mutants (Figure 6, B and C, right; P < 0.01). Transcriptome analysis revealed that nearly half of the HS-responsive genes were not induced or repressed in Mphsfa1U (Supplemental Figure S7, H and I). Some other HS-responsive genes were still induced in the mutant, but at a lower level (Supplemental Figure S7J; P < 0.001). These transcriptomic changes in Mphsfa1U were highly similar to those in Mphsfa1V (Figure 6, D and E). Consistent with Mphsfa1V (Figure 3C), transcriptomic defects in Mphsfa1U were related to the GO terms heat response and protein folding (Figure 6F) and oxidative stress response and fungal response (Supplemental Figure S7K). Likewise, HS-responsive metabolic changes in Mphsfa1U were comparable to those of Mphsfa1V (Figure 5B and Figure 6G). In general, the Mphsfa1U and Mphsfa1V mutants exhibited greater changes in the accumulation of various amino acids and sugars after 5-h of HS treatment compared to BC3 and Tak-1 (Figure 6H, Supplemental Figure S7L). These results indicate that MpHSFA1 plays comparable roles in the transcriptomic and metabolic regulation of HS responses.

Overexpression of MpHSFA1U or MpHSFA1V inhibits thallus branching and enhances thermotolerance

MpHSFA1 knockouts resulted in growth retardation, a possible tradeoff between growth and stress responses (Figure 3). We next examined if the overexpression of MpHSFA1U or MpHSFA1V (Supplemental Figure S8A) alone without HS would cause any developmental anomalies. The thalli of MpHSFA1 overexpression (ox) lines were four- to five-fold smaller than BC3 and Tak-1 thalli (Supplemental Figure S8B and S8C; P < 0.01). In addition, MpHSFA1 ox had longer plastochrons, causing the thallus to stretch, with fewer branches rather than appearing rounded (Supplemental Figure S8, C and D). Gender-based differences were observed in MpHSFA1-induced necrosis. Both MpHSFA1U and MpHSFA1V ox lines in the BC3 background exhibited necrosis of the growing points beginning on Day 14, while MpHSFA1 ox lines in the Tak-1 background showed no necrosis even on Day 28 (Supplemental Figure S8E). Despite the gender-based differences in necrotic phenotypes, all MpHSFA1 ox lines showed enhanced thermotolerance in both the BC3 and Tak-1 backgrounds (Supplemental Figure S8, F and G).

Transcriptomic profiles of MpHSFA1 ox lines under non-HS conditions would reflect the developmental and other biological processes controlled by HS. A comparison of such transcriptome data revealed roughly 100 genes commonly regulated by MpHSFA1U and MpHSFA1V overexpression and HS in both the BC3 and Tak-1 backgrounds (Supplemental Figure S9, A and B, Supplemental Data Set 2). However, their levels of induction and repression in response to MpHSFA1 ox were significantly smaller than those regulated by HS in the parental wild-type plants (Supplemental Figure S9, C and D). GO terms associated with temperature response and protein folding were consistently enriched in the overlapping DEGs (far right columns, Supplemental Figure S9E), although several HS-responsive genes were specifically induced in Tak-1 or BC3 plants (far left columns, Supplemental Figure S9E). Interestingly, HSEs were highly enriched in HS-responsive genes in both BC3 and Tak-1 (far left columns, Supplemental Figure S9F), but not in the genes controlled by MpHSFA1 ox under non-HS conditions. Taken together, MpHSFA1 overexpression alone only partially induced HS responses and transcriptomic changes and did not fully reproduce thermal responses, which might require the activation of additional thermal signaling pathways.

Gender-specific developmental and thermotolerance phenotypes associated with MpHSFA1U and MpHSFA1V

There was no clear gender-based differences in transcriptomic changes caused by HS (Supplemental Figure S7, A and B) or MpHSFA1 overexpression (Supplemental Figure S8) when we compared Tak-1 with BC3, which is a female introgression line of Tak-1. Nevertheless, BC3 plants were more sensitive to HS than Tak-1 plants (Figure 7, A and B; P < 0.001). To further explore the differences in thermal sensitivity, we compared the transcript levels of sex chromosome genes using fragments per kilobase of transcript per million mapped reads (FPKM) values as a measure of expression levels (Figure 7C and Supplemental Figure S10A). Three and 12 genes out of 75 U chromosome and 99 V chromosome genes showed altered transcript levels after HS treatment (FDR < 0.1), but this was independent of gender and Mphsfa1 knockout (Supplemental Figure S10A). Interestingly, however, the FPKM value of MpHSFA1V was 2- to 2.5-fold higher than that of MpHSFA1U under non-HS conditions (Figure 7C). In addition, genes downregulated specifically in Tak-1 plants, but not in BC3 or Mphsfa1 plants, were associated with DNA conformation changes (Supplemental Figure S10B), which may explain the gender-based differences in HS sensitivity (Ray et al., 2019; Montgomery et al., 2020).

Figure 7

Male and female Marchantia plants exhibit different levels of sensitivity to HS, and gender-specific reproductive organ development. A and B, Thermal tolerance of 7-day-old BC3 and Tak-1 plants after HS treatment at 37°C for 40 h. Representative plant images after 5 days of recovery are shown in A. Ratios of surviving to total plants after 37°C treatment for 24 and 40 h are plotted in boxplots (B) with the median values and the first and the third quartiles. ***P < 0.001 determined by two-tailed Student’s t test; n = 80 gemmalings of each genotype from four independent experiments. C, FPKM values for HSFA1V and HSFA1U expression under control conditions (Cntr), HS for 1 and 5 h in Tak-1 and BC3 plants, respectively. Color key represents log2(FPKM +1). D and E, Heat sensitivity of antheridiophores in Tak-1 and Mphsfa1V and Mphsfb1 mutants. Representative photographs in (D) were taken 15 days after transfer to far red light, followed by HS at 37°C for 24 h and then recovered at 22°C for 3 days. The white arrows indicate the damaged antheridiophores. The numbers of antheridiophores formed under far red light at 9, 11, and 13 days are shown in (E). The heatmap shows the mean values of 10 plants (n = 10) of each genotype. Color key represents the number of antheridiophores. ***P < 0.001, calculated by Mann–Whitney U test. F and G, Heat sensitivity of archegoniophores in BC3 and Mphsfa1U plants. Archegoniophore formation was induced by far red light for 24 days and HS at 37°C for h, followed by 3 days of recovery at 22°C. The white arrows indicate the damaged archegoniophores. G, Number of archegoniophores after 18, 20, and 22 days of growth under far red light. The mean numbers of archegoniophores observed at each time point are shown in the heatmap. Color key represents the number of archegoniophores. **P < 0.01, ***P < 0.001, calculated by Mann–Whitney U test. H, Summary of HSFA and HSFB evolutionary history and function.

Finally, we analyzed reproductive organ development as a haplotype-specific biological process. Both male and female Mphsfa1 mutants showed delayed antheridiophore or archegoniophore development compared to Tak-1 (Figure 7D (Cntr) and 7E) or BC3 (7F (Cntr) and 7G). Mphsfa1V and Mphsfa1U plants also exhibited more yellowish stems or reproductive organs after HS (Figure 7, D [HS] and 7, F [HS], white arrows), pointing to their roles in the development and thermotolerance of reproductive tissues. Surprisingly, Mphsfb1 plants produced no antheridiophores under the same condition (Figure 7, D and E). Based on this observation, together with the results of transcriptomic analysis (Figure 3), we propose that MpHSFA1 controls reproductive organ formation as a part of a HS response, whereas MpHSFB1 directly regulates antheridophore development. The data collectively suggest that MpHSFA1 is a master regulator of HS responses, while MpHSFB1 has more specialized functions in addition to temperature responses.

Land plants often face gradual temperature changes over a long period of time, which could also be regulated by MpHSFA1 and MpHSFB1. We therefore treated plants at 22°C, 29°C, and 32°C and monitored their growth over 21 days (Supplemental Figure S11, A–E). Mphsfa1V plants were hypersensitive to 29°C, as shown by the significantly reduced plant size (P = 0.011) at this temperature compared to 22°C (Supplemental Figure S11F). In contrast, the thallus growth of Tak-1 and Mphsfb1 was similar between 22°C and 29°C. At 32°C, thallus growth was arrested after 5 days of treatment regardless of genotype. These results suggest that MpHSFA1 is responsible for controlling wider temperature responses, spanning ambient temperature changes to HS responses, while MpHSFB1 plays important roles preferentially at higher temperatures.

Gemma cup formation was also promoted at 29°C in Tak-1, but not in Mphsfa1V, which had a similar number of gemma cups at both 22°C and 29°C (Supplemental Figure S11G; P = 0.018). Female BC3 plants showed a severe growth defect (P < 2 × 10−16) at 29°C, where their thalli tended to curl and burrow into the agar medium (Supplemental Figure S12, A–C). The growth defect was more severe in Mphsfa1U, which were smaller at 29°C (Supplemental Figure S11; P = 0.021). The observed gender-associated difference aligns with the lower survival rate observed in BC3 plants compared to Tak-1 plants (Figure 7B).

Discussion

Terrestrialization of land plants is generally viewed as the most important event in the natural history of life. Successful terrestrialization required drastic changes in the morphological features and intracellular mechanisms of plants to help them adapt to harsh terrestrial environmental conditions, such as high temperature and water deficiency, a stark contrast to the likely stable aquatic environment. The fixation of HSFs and the subsequent rapid evolution of the HSFB subclass occurred during the early evolution of Streptophyta plants (Figures 1 and 2). This evolutionary mode contrasts to that of typical TFs including HSFA subclass members within the vascular plant clade (Scharf et al., 2012), in which the diversification of TFs is gradual, as supported by phylogenetic trends with small phylogenetic branch lengths. The saltation of protein sequences appears as an unbalanced long branch length in a phylogenetic tree. Such long branches, which were described previously in plant TF phylogenies (D), could have also augmented the evolution of ancestral green species to widen the physiological and developmental functions that are inherited by extant land plant species. The model for plant HSF diversification and functional expansion is summarized in Figure 7H.

The rapidly evolved HSFB subclass control morphological processes under non-HS conditions (Figures 5, 7, D and E) and the expression of various genes related to development or cell wall biosynthesis (Figures 3 and 4). The developmental defects in Mphsfb1, such as defective antheridophore formation (Figure 7, D and E), could be mediated by auxin signaling pathways. The vegetative-to-reproductive phase transition in angiosperms requires auxin accumulation (Cheng et al., 2007; Krizek, 2011). Analysis of organ-specific transcriptomes in M. polymorpha suggested the co-expression of MpHSFB1 with auxin biosynthesis and perception genes such as TAA, YUCCA, and AUX/IAA in reproductive organs (Flores-Sandoval et al., 2018). These lines of evidence prompted the suggestion that auxin mediates some HSF-dependent developmental processes in M. polymorpha and in A. thaliana. Interestingly, genes for abiotic stress responses, defense responses, and acquired system resistance were anomalously expressed in Mphsfb1 even under non-HS conditions (Figure 3). Similar defense-related functions of HSFB were observed in S. lycopersicum, where an SlHsfb1-RNAi mutant showed altered proteome profiles associated with abiotic stress responses, defense responses, and cell-wall modification (Fragkostefanakis et al., 2019). In addition, Arabidopsis AtHSFB1 is known to regulate the growth-to-defense transition upon pathogen invasion (Pajerowska-Mukhtar et al., 2012). Hence, the defense-related functions of HSFB, as well as their developmental functions, were evolutionarily fixed prior to the separation of vascular and nonvascular plants.

Our transcriptome and metabolome analyses revealed the accumulation of amino acids such as proline, alanine, phenylalanine, and tyrosine under HS (Figure 4). Proline serves as a general stress protectant (Hare and Cress, 1997), while the aromatic amino acids phenylalanine and tyrosine are two key compounds used to produce phenylpropanoid-derived secondary metabolites (Vogt, 2010). The genes for phenylpropanoid pathways were enriched under HS in multiple species including Klebsormidium (labeled by the term cell structure in Supplemental Figure S5E). In addition, similar changes in metabolites and metabolic pathways were previously observed in heat-stressed Arabidopsis, tomato, and green algae (Fragkostefanakis et al., 2016; de Vries et al., 2020; Paupière et al., 2020). These observations indicate that the accumulation of phenylpropanoid compounds is a well conserved HS response inherited from aquatic green species. Genes involved in energy production were also enriched under HS in multiple species, indicating the important role of primary metabolic changes that were conserved from aquatic to terrestrial plants (Supplemental Figure S5E). The accumulation of components in both the phenylpropanoid and energy pathways (amino acids, sugars, and organic acids) is regulated by both MpHSFA1 and MpHSFB1 (Figure 4). While the metabolomic profiles of HSF-null mutants in angiosperms are yet to be explored, our findings suggest that the metabolic functions of HSF existed in ancestral green algae and were retained in the HSFA and HSFB subclasses after the split.

This study also focused on gender-based differences in HS responses, which were not explained by the sequence diversity of MpHSFA1U and MpHSFA1V (Figure 7). Instead, transcriptional regulation of MpHSFA1V and MpHSFA1U may partially account for the differences (Figure 7). For example, MpHSFA1V (male), with high steady-state expression levels, contains 80 TF-binding motifs, including 4 HSEs, whereas MpHSFA1U (female), which is expressed at low levels, only contains 41 motifs with no HSE (within the 2-kb promoter region). MpHSFA1U and MpHSFA1V share CDS- and protein-sequence identities of 63.5% and 74.5%, respectively, with different exon-intron structures. We also found that HSFA1U gene is located on the U chromosome in Marchantia inflexa, which separated from M. polymorpha more than 100 million years ago (Villarreal et al., 2016). However, since BC3 is an ∼90% introgression line, there might be few SNPs that contribute to the different HS responses of male and female M. polymorpha (Bowman et al., 2017) While this concept must be confirmed experimentally, the long-term retention of the two HSFA genes on Marchantia sex chromosomes further supports the regulatory or functional distinctions between the gender-specific HSFAs.

Collectively, this study focused on HSFA and HSFB in M. polymorpha and revealed their common and subclass-specific functions (Figures 3–7). While this and other previous studies mainly focused on the functions of HSF genes in HS responses, our transcriptome and cis-regulatory element analyses strongly suggested the interplay of these HSFs with different types of TFs under HS. For example, besides the cis-regulatory motifs (e.g. HSE) that could be found in commonly regulated genes in Mphsfa1 and Mphsfb1, cis-regulatory motifs for HSFs and stress-related TFs were highly enriched in Mphsfa1-specific genes, while WRKY-binding motifs were highly ranked among Mphsfb1-specific genes (Supplemental Figure S3 and S4). The different coupling partners would have given rise to functional distinctions between GRNs downstream of HSFA and HSFB either directly or indirectly. Duplicated genes are often maintained in the genome for a long period of time once one of the duplicates has acquired new important functions. This consequently leads to the expansion of GRNs (Voordeckers et al., 2015). The subsequent functional distinctions between the HSFA and HSFB subclasses would exemplify typical evolutionary processes that occurred in early plant life.

The GRNs underlying HS responses are initiated when land plants and algae face rapid and drastic temperature increases. These GRNs are centered around HSFs. A single canonical HSF (CrHSF1) in Chlamydomonas reinhardtii is a key regulator of HS that directly activates HSPs to protect cells from mis-folded proteins (Schmollinger et al., 2013). CrHSF1 possesses features typical of Class-A HSFs, suggesting their conserved roles across algae and land plants (Schroda et al., 2015). In addition, the expression of HS marker genes such as HSP70, HSFP90, and sHSPs was induced at 15 min of HS treatment and lasted for 2 h in C. reinhardtii and K. nitens (Monte et al., 2020), which is similar to our observations in A. thaliana seedlings or M. polymorpha gemmalings at 36°C to 42°C (Supplemental Figure S4; Larkindale and Vierling, 2008; Myouga et al., 2006). This indicates that rapid transcriptional responses upon HS are similar in green algae and young land plant seedlings. Considering that the aquatic habitat is subject to prolonged HS due to the high heat capacity of water, transcriptomic comparisons of prolonged HS responses between aquatic and land plants would be of interest in the future. MpHSFA1 and MpHSFB1 likely control plant responses at different temperature ranges (Supplemental Figures S11 and S12). The activation of HSFs and downstream TFs or genes at different temperature ranges may thus be a part of their functional distinctions.

The Streptophyta alga K. nitens and the bryophyte M. polymorpha contain two HSFs belonging to Class A and Class B. Our sequence analysis revealed that some charophyte algae lack HSFA or HSFB (Figure 2), implying that the two classes of HSFs were not evolutionarily fixed in some algal lineages. Our study showed that MpHSFA1 is a key regulator of HS response, similar to CrHSF1 and other Class A HSFs in vascular plants. On the other hand, MpHSFB1 possesses functions distinct from Class A HSFs as well as Class B HSFs in land plants. Given the divergent functions demonstrated in Marchantia, Arabidopsis, and tomato (Ikeda et al., 2011; Fragkostefanakis et al., 2019), HSFB1 likely acquired specialized functions during the evolution from algae to land plants to extend HS GRNs.

Materials and methods

Growth conditions and morphological characterization

Marchantia polymorpha wild type (accession BC3 [female] and Tak-1 [male]), MpHSFA1 overexpression lines, and Mphsfa1 and Mphsfb1 knockout mutants were grown on 1/2 × Gamborg B5 medium (pH 5.5) with 1% agar under continuous light (50–60 µmol m−2s−1, 6500K Cool Daylight LED, Philips) at 22°C. For morphological characterization, gemmae were grown under a 16-h light/8-h dark photoperiod for 10 days at 22°C, and the gemmalings were then transferred to continuous light at 22°C. Plant images were taken daily from Days 10 to 21 for the analysis of plant area, number of thalli, and aspect ratio of the ellipse covering the thallus (ratio of the length of the shorter axis to the longer axis). These morphological characters were extracted with OpenCV2 (Pulli et al., 2012).

HS treatment of M. polymorpha

Ten to 20 wild-type and mutant gemmae from the respective hsfa1 lines were grown on 1/2 × B5 (pH 5.5) medium with 1% agar at 22°C under a 16-h/8-h light/dark period for 10 days. The 10-day-old plants were treated at 37°C for 24-h with continuous light in a growth chamber and then recovered for 14 days at 22°C. Plants were imaged on Days 1 (right after treatment) and 14. Survival ratio was calculated as the green pixel area/total thallus area. The green and thallus pixels are defined as the H, S, and V values between [30, 20, 255] and [80, 255, 255], respectively. The proportion of gemmalings that remained green on Day 1 and Day 14 were recorded for each genotype. To analyze the tolerance to harsh temperature changes, hsfb1, HSFA1U and HSFA1V ox mutants grown in agar plates were immersed in a water bath preheated to 37°C. Seven-day-old plants grown on 1/2 × B5 (pH 5.5) plates with 1% agar were immersed in a water bath at 37°C for 2 h (without light) and immediately transferred to a 22°C growth chamber with continuous light for recovery. Photographs were taken and survival rates were analyzed on Day 5 after HS.

For long-term HS, gemma from Mphsfa1V, Mphsfb1, Tak-1, Mphsfa1U, and BC3 plants were grown under a 16-h light/8-h dark photoperiod for 7 days at 22°C. The gemmalings were then transferred to growth chambers under a 16-h light/8-h dark photoperiod at 22°C, 29°C, or 32°C for 21 days. Plant images were taken daily for the analysis of plant area and number of gemma cups. Plant size was extracted with OpenCV2 (Pulli et al., 2012).

HSF overexpression and knockout mutant lines in M. polymorpha

All cloning and transformation were performed as previously described (Wu et al., 2021). For overexpression mutant lines, protein-coding sequences of MpHSFA1U and MpHSFA1V were amplified from BC3 and Tak-1 plants, respectively. The DNA fragments were inserted into the binary vector pMpGWB336 that expresses an N-terminal tagRFP-tagged protein under the control of the MpEF1α promoter. The plasmids were transformed into M. polymorpha thalli via Agrobacterium-mediated transformation (Ishizaki et al., 2016), and successful transformants were screened on 1/2 × B5 medium plates with 100-µg mL−1 cefotaxime and 0.5-μM chlorosulfuron. Transgene expression was confirmed in T2 generation plants using qPCR, and two independent ox lines were selected based on the transcript level (Supplemental Figure S7A). Gene locus identifiers and primers used for sequencing and cloning are listed in Supplemental Data Set 5.

Mphsfa1U, Mphsfa1V, and Mphsfb1 knockout mutants were generated using the clustered regularly interspaced short palindromic repeats/CRISPR Associated protein 9 (CRISPR/Cas9) method. Two gRNAs flanking the short DNA fragment of the target gene were designed and cloned into the pMpGE011 vector. Agrobacterium-mediated transformation of thalli was performed as previously described (Kubota et al., 2013), followed by the screening of transformants with 100-µg mL−1 cefotaxime and 0.5-µM chlorosulfuron. For double transformation, thalli were co-transformed using two constructs with the same selectable marker. All experiments were conducted using the T2 generation. Detailed information about the Mphsfa1V, Mphsfb1, Mphsfa1Vb1, and Mphsfa1U knockout mutants and gRNA sequences is provided in Supplemental Figures S2G, S6E, and S7C.

RNA sequencing and data analysis

For M. polymorpha, gemmalings were grown in 1/2 × B5 liquid medium for 7 days under continuous light and subjected to HS at 37°C for 0, 1, and 5 h. Total RNA was isolated from the gemmalings using an RNeasy Plant Mini Kit (Qiagen, Hilden, Germany). For K. nitens, droplets of algal suspension from a fully grown plate were plated on C-medium plates and incubated for 1 month. The algal plates were then immersed into a water bath at 37°C for 2 h. Library preparation, including polyA enrichment and sequencing of 150-bp paired-end reads on the NovaSeq platform, was performed at Novogene, Singapore.

Raw sequences from M. polymorpha and K. nitens were mapped to the reference genomes (M. polymorpha genome v.5.1 and K. nitens GCA_000708835.1_ASM70883v1 from NCBI) using STAR (v.2.5.3a) (Dobin et al., 2013) with default settings. Differentially expressed genes (DEGs) were identified with the R package DESeq2 (Love et al., 2014), followed by comparisons to the expression level at 0 h with the threshold of |Log2-fold change| ≥ 1 and False Discovery Rate (FDR) < 0.1 using the “lfcShrink ashr” method. Arabidopsis orthologs of M. polymorpha and K. nitens genes were identified as previously described (Fischer et al., 2011) for GO analysis. GO term enrichment analysis was performed with the AgriGO v.2 database (bioinfo.cau.edu.cn/agriGO) using standard settings. Significantly enriched GO terms (FDR < 0.05) were visualized as dot plots using the R package ggplot2. The temperature used for RNA-sequencing was determined based on the induction levels of HSFA1V, HSFB1, and a heat shock protein gene (sHSP) at 22°C to 40°C (Supplemental Figure S3A). The three genes showed the highest expression at 37°C and decreased expression levels in the hsfa1V and hsfb1 mutants. For the analysis of sex chromosome genes, trimmed sequences from M. polymorpha were mapped to the M. polymorpha (v.3.1, marchantia.info) reference genome because the full gene annotations of both sex chromosomes were available in this version. The IDs were converted to v.6.1.

Motif enrichment analysis

AME from MEME suite (Bailey et al., 2009) was used to identify enriched TF binding motifs in the 2-kb region upstream from the transcription start site of genes assigned to each DEG group containing genotype-specific genes (Supplemental Figure S3D) or HS-induced TF genes (Figure 2C and Supplemental Figure S3E). The genes were compared with nonheat-responsive genes (Supplemental Figure S3D) or nonheat-responsive TFs (Figure 2C and Supplemental Figure S3E) to calculate the enrichment (adjust P-value) by Wilcoxon Rank sum test, adjusting by Benjamini–Hochberg (BH) post-hoc test. The parameters were set as follows: ‐‐verbose 1 ‐‐scoring avg ‐‐method fisher ‐‐hit‐lo‐fraction 0.25 ‐‐evalue‐report‐threshold 10.0 ‐‐control ‐‐shuffle‐‐ ‐‐kmer 2. The Arabidopsis DAP‐Seq database (O’Malley et al., 2016) was used as a reference.

RNA extraction, cDNA synthesis, and qRT-PCR

Seven-day-old gemmalings grown in liquid 1/2 × B5 medium were treated at 37°C for 0, 1, and 5 h with continuous light. After treatment, plant samples were harvested, immediately frozen in liquid N2, and stored at −80°C until further use. Total RNA was extracted from M. polymorpha gemmalings using an RNeasy Plant Mini Kit (Qiagen, Hilden, Germany). First-strand cDNA was synthesized using the GoScript Reverse Transcription System (Promega, USA). qRT-PCR was performed as previously described (Wu et al., 2020). Data were analyzed according to the manufacturer’s instructions using CFX Manager v3.1 (Bio-Rad, USA). Gene expression levels were normalized using the expression of M. polymorpha ACTIN 7 (MpACT7). The primers used are listed in Supplemental Data Set 5.

ChIP and ChIP-sequencing

Approximately 1.5 g of 7-day-old gemmaling tissue from Tak-1 and the MpHSFA1V overexpression lines were cross-linked with 1% formaldehyde in a vacuum chamber twice for 5 min. Subsequently, the chromatin DNA was isolated and sonicated as described previously (Wu et al., 2018). The DNA was immunoprecipitated with 25-µL anti-tagRFP magnetic beads (rtd-20, Chromotek, Germany) into 1 mL of solution containing 1.1% Trion X-100, 1.2-mM EDTA, 16.7-mM Tris–HCl, pH 8.0, 167-mM NaCl, and protease inhibitor (cOmplete Protease Inhibitor Cocktail, Roche) (Wu et al., 2018). The precipitated DNA was purified with a Qiaquick PCR purification kit (Qiagen, Germany) and processed further for ChIP-sequencing. Library preparation and sequencing of 150-bp paired-end reads were performed on the NovaSeq platform in Novogene, Singapore. Two biological replicates of Tak-1 and MpHSFA1V overexpression lines were included. Bowtie2 (v2.3.5) software (Langmead and Salzberg, 2012) was used for mapping using default settings, except that we discarded reads matching multiple loci, as they might introduce error signals by repeat counting. Peak calling was performed with the model-based analysis software MACS2 (Zhang et al., 2008) using Tak-1 data as a control. Two biological replicates were peak-called separately, and only the peaks that overlapped in both replicates with an irreproducible discovery rate <0.05 were included in further analysis. The ChIPseeker (Yu et al., 2015) package in R was used for downstream analyses including peak annotation and genomic feature visualization. GO analysis was carried out using the AgriGO v.2 database. Motif discovery analysis was performed in AME from MEME suite (Bailey et al., 2009). The enriched motifs were compared with the DAP-seq database (O’Malley et al., 2016).

Interspecies comparison of HS-induced transcriptomic changes

A publicly available microarray dataset was retrieved from GEO (Supplemental Figure S1A; Liu et al., 2011; Liu and Charng, 2013; Zhang et al., 2013; Mishra et al., 2016; Cortijo et al., 2017; Keller et al., 2017; Albihlal et al., 2018; Nuño-Cabanes et al., 2020). The microarray data were processed with quantile normalization and background correction using the robust multi-array average (RMA) with the R packages “affy” and “gcrma” (Wu et al., 2004; Gautier et al., 2004). The “limma” package (Smyth, 2005) was used to compute log2-transformed values of transcripts in a linear model. DEGs from the corresponding control group were identified using log2-fold change (FC) value >1 and FDR <0.05 as the threshold. Affymetrix probe IDs were then converted to AGI locus identification numbers on the BAR website (Waese et al., 2017). Probe IDs with no match or ambiguously matching multiple loci were discarded.

RNA-seq datasets for S. cerevisiae, A. thaliana, O. sativa, and S. lycopersicum were obtained from the GEO database (Supplemental Figure S1A). The raw fastq files were mapped to the respective genomes using STAR mapper to generate raw read counts. The datasets were first filtered by array quality weight. Datasets with moderated-t-statistics that were out of range of other datasets were eliminated. All data were processed using the package LIMMA (Smyth, 2005) to minimize the variations caused by different algorithms used for data pre-processing. Affymetrix and Agilent microarray data were background corrected using the Robust MultiChip Average (RMA) method, while RNA-Seq data were treated with voom, where the mean–variance relationship of the log-read counts was estimated and given a precision weight (Law et al., 2014). The data were then normalized using quantile normalization and converted into log-scale before entering the limma empirical Bayes analysis pipeline. To access the data quality after processing, the distance between array, array intensity distribution, and mean–variance trend were evaluated.

The genome versions used in this study were as follows: S. cerevisiae S288C, A. thaliana-TAIR10, K. nitens-GCA_000708835.1_ASM70883v1, S. lycopersicum-SL 3.0, M. polymorpha-v3.1 and v5.1, O. sativa-IRGSP 1.0.

Orthologous gene analysis

Whole proteome sequences of S. cerevisiae, K. nitens, A. thaliana, O. sativa, S. lycopersicum and M. polymorpha were downloaded in FASTA format from the following databases: S. cerevisiae: S288C from Saccharomyces genome database (https://www.yeastgenome.org/), K. nitens: GCA_000708835.1_ASM70883v1 from NCBI (https://www.ncbi.nlm.nih.gov/), A. thaliana: TAIR10 from TAIR database (https://www.arabidopsis.org/), O. sativa and S. lycopersicum: IRGSP 1.0 and SL3.0 from Phytozome v13 (https://phytozome-next.jgi.doe.gov/), and M. polymorpha: v.5.1 from MarpolBase (https://marchantia.info/).

The proteome sequences from primary transcripts were used to construct OrthoMCL groups by sorting the sequences to pre-determined orthogroups (v6.1) with BLASTp (https://orthomcl.org/orthomcl/) (Fischer et al., 2011). The full list of orthogroups in each species is provided in Supplemental Data Set 1. In total, we identified 4,767 (S. cerevisiae), 11,246 (K. nitens), 7,804 (M. polymorpha), 24,794 (O. sativa), 14,689 (A. thaliana) and 10,892 (S. lycopersicum) orthogroups in six species from a total of 495,339 pre-defined orthogroups in OrthoMCL v6.1. The genes assigned to the OrthoMCL groups were used for further analysis. The HS responses of these orthologous genes were analyzed and presented as the degree of overlap (calculated as |A∩ B|/min(|A|, |B|), where |A| and |B| represent the number of stress-related genes in stress A and B, respectively; Ma et al., 2014) of orthologous genes between species.

Phylogenetic analysis

HSF DBD domain sequences from A. thaliana and Oryza sativa spp. Japonica, Picea abies, Ginkgo biloba, Azolla filiculoides, Salvinia cucullate, Selaginella moellendorffii, Ceratodon purpureus (transcriptome), Physcomitrella patens, M. polymorpha, Spirogyra pratensis, Coleochaete orbicularis, Chara braunii, Mesotaenium endliicherianum, and Klebsormidium flaccidum were collected from the Tapscan database and JGI PhycoCosm (https://phycocosm.jgi.doe.gov/phycocosm/home) (Wilhelmsson et al., 2017). HSF sequences from Chlorophyta and Rhodophyta species; Chlamydomonas reinhardtii, Volvox carteri, Chlorella sp. NC64A, Chlorella vulgaris, Porphyra umbilicalis, Cyanidioschyzon merolae, Chondrus crispus, and Galdieria sulphuraria were included as outliers to construct a phylogenetic tree. The HSF sequences were preprocessed by removing redundant sequences that showed 90% or higher identity to other sequences. The remaining sequences were aligned using MAFFT version 7, followed by the removal of gapped positions containing ≥10% gaps. MpHSFA1U, MpHSFA1V, and MpHSFB1 sequences were manually added to the multiple sequence alignment, as they were removed due to high similarity to other entries. The nonredundant set of HSF sequences (110 sequences with 89 positions, Supplemental Files 1 and 2) was used for phylogenetic analysis by IQ-TREE multicore version v2.1.2. The phylogenetic tree was constructed with the LG + R5 substitution model and 1,000 ultrafast bootstrap iterations. The model was selected based on Bayesian Information Criterion. For duplication length analysis (Vosseberg et al., 2021), phylogenetic distance to the duplication node uniting HSFA and HSFB branches was calculated for each HSFA or HSFB sequence. This distance was divided by the stem length of HSFA or HSFB branches (0.083 or 0.1924, respectively, Figure 2A). The duplication length was negatively log-transformed and presented in a bar histogram with frequency density. The common ancestor sequence was reconstructed using phyml (Guindon et al., 2005).

Structure homology modeling

The structures of MpHSFA1 and MpHSFB1 DBD were modeled using Modeller software version 9.23 (Eswar et al., 2006). A template atomic structure of animal HSF1 was obtained from the Protein Data Bank (PDB, ID = 5d5u; Neudegger et al., 2016). The protein sequences of MpHSFA1 and MpHSFB1 DBD were aligned with the HsHSF1 sequence in Modeller. Five structures each of HSFA1 and HSFB1 were modeled automatically with default settings. Modeled structures were assessed by the discrete optimized protein energy (DOPE) score, and the structures with the lowest DOPE values were used for further analyses. The PyMOL program was used to color atoms and generate the final images.

Entropy measurement at each position in the MSA

The Shannon entropy and JSD scores were used to estimate the conservation of residues for each position in plant HSF sequences (108 sequences with 89 positions, see above for multiple sequence alignment). Shannon entropy and JSD scores were calculated with the BLOSUM62 substitution matrix for estimating substitution probabilities and the Protein Residue Conservation Prediction program for background frequency of residue variability (Capra and Singh, 2007). The JSD score is presented in pseudo-color scale on the HSF structures. The Shannon entropy at each position is presented in Seq-Logo format.

Metabolite extraction and metabolic profiling by GC–MS

Seven-day-old M. polymorpha gemmalings (∼50 mg) were collected and immediately frozen in liquid nitrogen prior to storage at −80°C right after HS treatment at 37°C for 0, 1, and 5 h. The samples were then homogenized using a mortar and pestle precooled with liquid nitrogen and extracted in 700-µL methanol, and 60-µL of internal standard (0.2-mg ribitol mL−1 H2O) was subsequently added as a quantification standard. The extraction, derivatization, standard addition, and sample injection were conducted as described previously (Lisec et al., 2006). Untargeted GC–MS analysis was performed using a previously described system (Avin-Wittenberg et al., 2015). In brief, samples were separated using a DB-35ms Ultra Inert column (30 m × 0.25 mm × 0.25 µm, Agilent, USA) with Agilent 7200 Q-TOF GC–MS coupled to Agilent 7890B Gas Chromatograph by Chemical, Molecular, and Materials Analysis Centre (CMMAC), National University of Singapore. Raw MS data were extracted using the MassHunter Profinder suite (B.08.00, Agilent), and resulting peaks were aligned and mapped to the NIST database (2017) using Mass Profiler Professional software (v 8.0, Agilent) with default setting. Student’s t test (P < 0.05) was applied to generate differentially accumulated compounds in HS versus control conditions in at least one pair-wise comparison.

Statistical analysis

Statistical analysis and graph construction were performed in R (v 4.0.4). The results of statistical analysis are provided in Supplemental Data Set 6.

Accession numbers

Sequence information from this article can be found in the Marchantia database (https://marchantia.info/) or under the following accession numbers: MpHSFA1U: MpUg00290, MpHSFA1V: MpVg00470, and MpHSFB1: Mp4g12230.

Raw sequencing data are available at the NCBI website under accession number GSE178776.

Supplemental data

The following materials are available in the online version of this article.

Orthology analysis of HS-responsive genes.

Structural analysis of HSFA1 and HSFB1 and identification of hsfa1 and hsfb1 mutants in Marchantia.

Transcriptomic changes and putative cis-regulatory elements (pCREs) in response to HS in Marchantia.

Analysis of HSE-containing genes in Mphsfa1V and Mphsfb1 mutants.

Binding motifs of MpHSFA1.

Characterization of Mphsfa1V and Mphsfb1 mutants.

Identification and characterization of Mphsfa1U mutants.

Characterization of HSFA1U and HSFA1V ox lines in Marchantia.

RNA-seq analysis of HSFA1ox lines in the Tak-1 and BC3 backgrounds.

Transcriptomic changes in genes in BC3 and Tak-1 plants and in U and V chromosomes.

Sensitivity of Mphsfa1V and Mphsfb1 mutants at 29°C.

Sensitivity of Mphsfa1U mutants and BC3 plants at 29°C.

Orthogroups and orthologous HS-responsive DEGs from the selected taxa.

Differentially expressed genes in wild type and hsf mutant plants.

Metabolite contents under HS and control conditions.

Genes differentially expressed between three wild-type accessions of M. polymorpha (Tak-1, Tak-2, and BC3).

Primers used in this study.

Results of statistical analysis.

Sequence alignment of HSF proteins used for phylogenetic analysis.

Phylogenetic tree file for Figure 2A.

Funding

This study was supported by the Agency for Science, Technology and Research (A*STAR) Singapore under the industry alignment fund pre-positioning program; High Performance Precision Agriculture (HiPPA) system (A19E4a0101) and by the Singapore-MIT Alliance for Research and Technology; Disruptive & Sustainable Technologies for Agricultural Precision (DiSTAP). We thank Drs Nahiko Ohama (TLL) and Ming-Jung Liu (ABRCST, Taiwan) for their valuable comments and discussion.

Conflict of interest statement. The authors have no conflicts of interest.

Click here for additional data file.