Transcriptional acclimation to warming temperatures of the Australian alpine herb Wahlenbergia ceracea

Understanding the molecular basis of plant heat tolerance helps to predict the consequences of a warming climate on plant performance, particularly in vulnerable environments. Our current understanding comes primarily from studies in Arabidopsis thaliana and selected crops exposed to short and intense heat stress. In this study we sought to characterise the molecular responses of an Australian alpine herb (Wahlenbergia ceracea) to growth under sustained moderate warming. We compared responses of pre-defined tolerant and sensitive lines, based on measures of photosynthetic stability, to growth under cool or warm temperatures to identify the pathways involved in thermal tolerance and acclimation. Under warmer growth temperatures, W. ceracea up-regulated genes involved in RNA metabolism, while down-regulating those involved in photosynthesis and pigment metabolism. In tolerant lines, genes related to photosystem II, light-dependent photosynthetic reactions, and chlorophyll metabolism were more strongly down-regulated. This suggests that the regulation of electron transport and its components may be involved in thermal acclimation. Our results also highlight the importance of hormonal gene networks, particularly ethylene, during longer-term moderate warming. In conclusion, our results point to post-transcriptional processes and the stabilisation of the electron transport chain as candidate mechanisms for thermal acclimation in W. ceracea. The study also revealed many non-orthologous temperature-responsive genes, whose characterization may enhance our understanding of physiological acclimation and have relevance for the biotechnological improvement of threatened species and crops.


Introduction
Extreme environments, such as alpine regions, are among those predicted to be most severely impacted by global warming (Gobiet et al., 2014), but native and wild plant species from these environments remain under-represented in genomic studies (Geange et al., 2021).
Transcriptomics and bioinformatics have enormous potential to reveal the genomic architecture of traits under selection, their variation, and plasticity and adaptation to climate change (Oomen and Hutchings, 2017;Oomen and Hutchings, 2022); to identify the genetic basis of physiological attributes and acclimatory responses that may underlie resilience to harsh climates (Rashid et al., 2020); using transcriptomic profiles is likely to yield fundamental advances and potential application for conservation planning (Hansen, 2010;Oomen and Hutchings, 2022).Therefore, application of these techniques to more species and natural populations is warranted.
Plant molecular responses to heat stress have been extensively explored in model plants and crops, such as Arabidopsis thaliana and rice (e.g., Niu and Xiang, 2018;Zhang et al., 2021;Li et al., 2021).Transcriptomic studies in Arabidopsis and other species revealed that responses to heat involve the regulation of genes encoding proteins collectively known as heat-shock proteins (HSPs), predominantly chaperones and transcription factors (Richter et al., 2010).The molecular chaperones, such as HSP70, HSP90, HSP100, and small HSPs (sHSPs), bind unfolded proteins preventing their aggregation and, in some cases, induce their correct re-folding or degradation (Boston et al., 1996;Richter et al., 2010;Wang et al., 2004).The heat-shock transcription factors (HSFs) regulate the expression of the HSPs in response to many abiotic stresses, including heat (Andrási et al., 2021;Haider et al., 2021a;Haider et al., 2021b).In addition, many phytohormones, such as abscisic acid (ABA) and ethylene, influence the responses to high temperatures (Li et al., 2021).The combination of such molecular responses helps maintain cellular function thereby promoting heat tolerance (Nievola et al., 2017;de Pinto et al., 2015;Yamori et al., 2014).
Most of these studies focus on short, intense heat events, often referred to as heat-shock (Wang et al., 2020), which can be used to infer the pathways that may underlie thermal acclimation.
Heat stress causes several changes to the physical and chemical properties inside cells, which in turn affect cellular function.Many proteins and enzymes lose their native folding thus impairing their catalytic function, with decreasing enzymatic activity even at temperatures slightly above their optimum for correct functioning (Daniel et al., 2008;Hobbs et al., 2013); these denatured proteins can form aggregates that eventually kill cells (Goraya et al., 2017;Lippmann et al., 2019).The concentration of reactive oxygen species (ROS) increases during heat stress, causing oxidative damage to proteins, DNA, and lipids (Suzuki and Mittler, 2006).High temperatures also increase membrane fluidity, altering their permeability to water, ions, and solutes, and affecting transmembrane proteins (Niu and Xiang, 2018;Cano-Ramirez et al., 2021).These changes in membrane stability can directly affect the function of photosystems and photosynthetic reaction centres, located on thylakoid membranes (Hugly et al., 1989).
Experimental heat-shocks may be useful as representations of heatwave events; however, plants also make metabolic and physiological adjustments to optimize their photosynthetic performance under more moderate or prolonged warming (Yamori et al., 2014;Way and Yamori, 2014;Larkindale et al., 2005).Responses to longer-term warming may involve distinct pathways compared to acute responses (Wang et al., 2018;Wang et al., 2020).For example, thermal acclimation to moderate or gradual warming in plants can be achieved by increasing the thermal optimum for photosynthesis, the electron transport rate, and/or stabilising Rubisco activase (Yamori et al., 2014;Way and Yamori, 2014).Thus, it is also vital that we study the molecular responses to longer-term moderate warming that plants must cope with through gradual, sustained rises in average temperature and that we examine these responses in non-model systems.
In this study, we provide de novo genome and transcriptome assemblies of the nonmodel alpine herb Wahlenbergia ceracea; we infer orthology with Arabidopsis thaliana for thousands of genes, thus providing preliminary annotations; and we investigate differential expression during thermal acclimation to moderate warming, summarising results via enrichment tests and prediction of transcription factors.This species was chosen for study since it demonstrates the capacity to acclimate to warmer growth temperatures, evident through measurable changes in photosynthetic stability (Notarnicola et al., 2021;Arnold et al., 2022), and because we had identified lines that differed significantly in their inherent thermal tolerance.Here, we conducted RNA-sequencing on leaf tissue of sensitive and heat tolerant lines to assess temperature-responsive gene regulation between tolerance categories.We see distinct signatures of regulation of posttranscriptional processes, hormonal (particularly ethylene) signalling, and stabilisation of electron transport chain during thermal acclimation when the sensitive and tolerant lines of this alpine plant are compared.

Generation of individuals for mRNA sequencing
Seeds of W. ceracea were collected in the field from mature capsules (F0), across the elevation range of the species (1590 -2095 m a.s.l.) from Charlotte's Pass to Mount Kosciuszko (NSW, Australia) and subjected to controlled crosses to generate F3 lines (Figure S1).A complete pedigree is presented in Table S1.
To identify tolerant and sensitive lines, we conducted a preliminary experiment where a total of 296 F3 individuals from 37 lines were grown in the glasshouse under cool (24/18C day/night) or warm (30/25C day/night) growth temperatures.Heat tolerance was assessed by measuring the temperature-dependent change in chlorophyll fluorescence (T-F0) as described in Arnold et al. (2021).We ranked lines using their average Tmax and selected the six F3 lines with the highest and lowest thermal tolerance respectively for the following experiment (tolerant lines average Tmax = 56.07°C;sensitive lines average Tmax = 50.83°C;Tmax [Tmax/tolerant -Tmax/sensitive] = 5.24°C; Table S2).
We randomly arranged three replicate plants for each line in three blocks (one replicate per block) that were equivalent in each of 2 capsules.Plants were bottom watered as needed and supplemented with liquid fertiliser every two weeks (Yates Thrive, Yates Australia).
Two weeks after transplantation, plants were subject to either cool (24/15C day/night, optimal for growth and flowering of post-seedling plants; Arnold et al., 2022) or warm temperatures (30/20C day/night).Plants grown at 30C have reduced fitness (seed number and viability) but are not sterile nor undergo photoinhibition [Fv/Fm (maximum PSII efficiency) > 0.70; Notarnicola et al., 2021].Experimental design: 12 lines × 3 replicates × 2 growth temperatures = 72 individuals (Figure S1).We let the plants grow for a further two weeks by which time they had produced new leaves acclimated to the respective temperature conditions.We took two ~1 mg leaf tissue per plant from the two youngest, fully expanded, healthy leaves.Leaves were harvested in a randomized order between 11:00 and 16:30, and immediately flash-frozen in liquid nitrogen and stored at -80C.Time of harvesting and developmental stage were recorded for each sample ('nonflowering' n = 56; 'early bud' n = 16).

Genome annotation and de novo transcriptome assembly
DNA was extracted from one F3 individual grown in the glasshouse.High-molecular weight DNA was extracted following Jones et al. (2020) and long-read sequencing was conducted on a PacBio Sequel II platform.We used Canu HiFi to assemble the PacBio reads, using a genome size estimate of 0.9 Gb and PacBio HiFi settings (Koren et al. 2017;Nurk et al. 2020).Reads mapping to non-plant contigs were identified using BlobTools and subsequently removed before assembling the remaining reads.Repeats were annotated and masked using RepeatModeler and RepeatMasker (Flynn et al., 2020).Gene annotation was conducted with the BRAKER2 pipeline with flag --epmode (Brůna et al. 2021).Gene predictions were guided by all eudicot orthogroups available in OrthoDB v.10 ( Kriventseva et al. 2019) and the mRNA sequencing dataset generated herein.For the latter, raw reads across all mRNA-seq samples were merged and aligned to the genome using Subjunc (Liao et al., 2013).The aligned reads were then included in the BRAKER2 predictions.We tested for completeness of the genome assembly and the annotations using Benchmarking Universal Single-Copy Orthologs (BUSCO; Seppey et al. 2019).We refer to the predicted gene annotations as the transcriptome.

mRNA-sequencing and analysis
Total RNA was extracted from snap-frozen and ground leaf tissue using the Spectrum TM Plant Total RNA Kit (Merck, Kenilworth, New Jersey, USA) with on-column DNase I digestion.
Purified RNA was quality checked for DNA contamination and integrity using semi-denaturing agarose-gel electrophoresis and the LabChip GX Touch (PerkinElmer).RNA was quantified with ND-1000 Spectrophotometer (NanoDrop Technologies) and Qubit 4 Fluorometer (RNA HS kit, Thermo Fisher Scientific).
Library preparation and sequencing were conducted by Deakin Genomics Centre.The construction of cDNA libraries was conducted with 250 ng RNA input using the Illumina stranded mRNA prep ligation protocol, as per manufacturer's instructions except that reaction volumes were scaled by one-half and 12 PCR cycles were used for library amplification (Illumina, San Diego, California, USA).Paired-end sequencing (2 × 150 bp) was conducted on the NovaSeq 6000 (Illumina).
Quantification was performed on transcripts defined from the de novo transcriptome assembly.
Transcript abundances were imported into R v.4.0.2 (R Core Team, 2021) and summarised to genelevel counts with Tximport, using the 'length scaled transcripts per million (TPM)' equation (Soneson et al., 2015).Only genes detected at > 1 count per million (CPM) in at least 18 samples were retained for analysis.We calculated scaling factors to normalise gene counts based on effective library size using the calcNormFactors function in edgeR (Robinson et al., 2010) using the Trimmed Mean of M-values (TMM) method.Differentially expressed genes (DEGs) were identified using the R package 'ThreeDRNAseq' (Guo et al., 2021) and 'limma' algorithms (Smyth, 2004).DEGs were determined using the limma.pipelinefunction in 'ThreeDRNAseq' with Benjamini-Hochberg P-value adjustment (adjusted P-value ≤ 0.05).Initial models included a factor combining the tolerance and temperature factors ('Tolerant_cool', 'Sensitive_cool', 'Tolerant_warm', and 'Sensitive_warm'; n per group = 18), as well as the developmental stage ('nonflowering' and 'early bud'), and collection time.However, preliminary analyses and Principal Component Analysis (PCA) did not reveal clustering of samples by developmental stage nor harvest time and so these factors were not included in the final models.Lines were treated as random effects: within-line correlation was estimated with the duplicateCorrelation (Smyth et al., 2005) function in 'limma'.Observational-level and sample-specific quality weights were assigned with the voomWithQualityWeights function (Liu et al., 2015).The full model was: gene expression ~ 0 + (tolerance_temperature), random effect = line.The exclusion of the intercept from the model and use of a combined factor allowed us to calculate differences in gene expression between predetermined contrasts (Table 1).
Enrichment tests were performed using Metascape.We conducted multi-gene list analyses separately for up-and down-regulated genes to compare enriched terms between tolerant and sensitive lines.These multi-gene lists included orthologs of significant genes in both tolerant and sensitive lines (Table S3, S4).We conducted gene annotation using default terms catalogues, except DrugBank was substituted with KEGG Pathway.We ran enrichment analyses with default settings (catalogues: GO Biological Process, KEGG Pathway, and WikiPathway; minimum overlap = 3; Pvalue cut-off = 0.01; minimum enrichment factor = 1.5).Transcription factors regulating the DEGs were predicted with TF2Network (Kulkarni et al., 2018) on up-and down-regulated genes separately, for tolerant and sensitive lines.

Genome and transcriptome assembly for Wahlenbergia ceracea
The size of the haploid genome assembled with Canu was 802.965Mb and accounted for 89% of the estimated genome size (900 Mb).The results for gene space completeness were comparable to other recent, high-quality plant genome assemblies; the assembly reported 91.5% of Embryophyte BUSCOs as complete, and 79.8% of BUSCOs as complete and single copy.
However, the assembly recorded low contiguity as the number of contigs was 1876, the contig N50 was 1.125 Mb, and the contig NG50 (based on the estimated genome size) was 0.985.The fragmentation of the genome assembly was related to the presence of repetitive sequences in transposable elements identified at the end of existing contigs.Gene annotations (transcriptome) with BRAKER2 predicted 223,827 exons for 87,730 proteins in the complete assembly.BUSCO testing revealed that the transcriptome was overall very similar to the genome assembly, although the former reported a marginal increase in fragmented BUSCOs and a corresponding decrease in complete single copy and missing BUSCOs.
Pseudo-alignment of mRNA-seq reads using Kallisto against the transcripts from the BRAKER2 prediction recovered 41,775 total transcript isoforms, which were summarised to genelevel counts with Tximport (Soneson et al., 2015) for a total of 37,307 genes.Of these, 22,585 were considered as expressed (≥ 1 CPM in ≥ 18 samples).On average, 69.7% of reads were pseudoaligned to the transcriptome per sample (Table S5).
Orthology inference revealed that W. ceracea shared the highest number of orthogroups with Manihot esculenta (9,357), Glycine max (9,275) and Vitis vinifera (9,239), among the species used for orthology inference (Table S6).We found 8,854 orthogroups between W. ceracea and Arabidopsis thaliana, and a total of 11,205 orthologs between the two species (Table S6, S7).
Overall, 81.4% ( 0.25 SD) of DEGs (see below) were orthologs with A. thaliana genes.Principal Component Analysis (PCA) revealed separation along PC2 (9.26%) reflecting the temperature treatment received (Figure 1a).In contrast, PC1 explained 14.69% of the variation but was not associated with any measured variable and may reflect genetic differences.To ensure we could observe a robust temperature response, we first verified that all W. ceracea lines were responsive to warmer conditions (contrast i).Indeed, we could detect 7,407 DEGs in response to temperature (Table 1, Figure 1b).Of these, 1,228 were detected when considering all lines (contrast i) but were not significant when tested within tolerant or sensitive lines (contrasts ii and iii) (Figure 1b), presumably due to reduced statistical power with smaller sample sizes.Overall, the results of differential expression and functional enrichment tests in contrast (i), presented in Figure S2 and Table S8, recapitulate those in contrasts (ii) and (iii), therefore we will only report and discuss results for the latter two contrasts.
A total of 4,732 genes (up = 2,305, down = 2,427) and 4,938 genes (up = 2,375, down = 2,563) were differentially expressed in response to temperature in tolerant and sensitive lines, respectively (Table 1).While 2,790 DEGs were in common between tolerance groups, 1,941 and 2,147 were temperature-responsive only in either tolerant or sensitive lines, respectively (Figure 1b).Despite the larger number of DEGs in sensitive lines, responses were more robust in tolerant lines, as shown by the smaller P-values (Figure 1c-d).Fold changes in expression are slightly skewed towards the left, suggesting higher changes in expression for downregulated genes.
In tolerant lines, the most significant upregulated gene in response to temperature did not have orthologs with A. thaliana and therefore may represent a novel gene; whereas the second and third genes were orthologs of bZIP28 (AT3G10800; Table S9), a transcription factor upregulated in response to heat that responds to unfolded proteins in the endoplasmic reticulum (Zhang et al., 2017).The most significantly downregulated DEG was ortholog of AT5G58490 (Table S9), described as a NAD(P)-binding Rossmann-fold superfamily protein.The second and third downregulated DEGs were orthologs of Phosphate Transporter 4;1 (PHT4;1, AT2G29650) localized to the thylakoid membrane and annotated as involved in transmembrane transport.
In sensitive lines, the most significant upregulated gene was an ortholog of AT4G35785 (Table S10), an RNA binding protein annotated as a positive regulator of mRNA splicing.The second and third top upregulated DEGs were respectively orthologous to AT4G36420, described as a ribosomal protein L12, and to NAP1-related protein (NRP2, AT1G18800), a protein that binds to histones and is involved in DNA recombination.The most downregulated gene was an ortholog of AT5G58490 (Table S7), also found downregulated in tolerant lines.The second and third most significant DEGs were respectively orthologous to AT2G21960, encoding a transmembrane protein involved in chlorophyll metabolic process and proteolysis, and to Vitamin E Deficient 1 (VTE1; AT4G32770), a cyclase involved in vitamin E synthesis.
To explore the function of temperature-responsive genes, we utilized functional information from orthologous genes in A. thaliana determined with OrthoFinder (Emms and Kelly, 2019).A total of 1,790 (upregulated) and 1,975 (downregulated) temperature-responsive genes in tolerant and sensitive lines (combined gene list from contrasts ii and iii; Tables S6, S7) were recognised and analyzed using Metascape.Overall, up-and down-regulated genes related to posttranscriptional regulation and photosynthesis, respectively (Figure 2a-b).Strikingly, the sensitive lines demonstrated a more significant enrichment of terms related to post-transcriptional regulation of gene expression among the up-regulated genes, such as ribonucleoprotein complex subunit biogenesis and RNA 3'-end processing.The opposite was observed for down-regulated genes, where tolerant lines exhibited a stronger enrichment for various terms including photosynthesis, thylakoid membrane organization, and NAD(P)H dehydrogenase complex assembly.
We next explored expression levels of the photosynthesis-associated DEGs to highlight those genes driving enrichment of this term and potentially involved in the higher Tmax values observed in tolerant lines.The GO and KEGG term 'photosynthesis' comprised 139 W. ceracea DEGs that were orthologs with 89 A. thaliana genes.These genes comprised many components of antennae and complexes involved in electron transport, such as two Protein Phosphatase 1 (PPH1), Evolving Complex 33 (OEC33), and ATP synthase C1 (ATPC1).All these genes were significant and more down-regulated in tolerant lines, while some of the antennae components were not significant in sensitive lines (Figure 3).Responses across all lines were consistent with group averages (Figure S3).
The 'protein folding' pathway, which includes HSPs, was significant among up-regulated genes only in contrast (i), but not in contrasts (ii) and (iii), suggesting no differential regulation between tolerant and sensitive lines.We found nine HSPs (five HSP90s, one HSP60, and three smallHSPs) and all were up-regulated under warmer temperatures (Figure S4, S5).

Gene regulatory networks regulating temperatureresponsive gene expression in heat tolerant and sensitive W. ceracea
Since we observed contrasting regulation between heat tolerant and sensitive W. ceracea, we interrogated whether these differences may be the result of transcription factor (TF) activity.
Using TF2Network to construct gene regulatory networks, we predicted 47 and 32 TFs from upregulated temperature-responsive DEGs in tolerant and sensitive lines, respectively, of which 26 were in common (Table S11).Many of the shared TFs predicted for both groups are involved in hormone signalling, particularly with ethylene (numerous Ethylene-Responsive Factors) but also chitin, jasmonic acid, cytokinin, and gibberellin.The predicted TFs that were not shared between tolerant and sensitive lines were also involved in ethylene signalling, suggesting a role of this hormone during temperature responses in W. ceracea.
Many more TFs were predicted for down-regulated genes; 190 in tolerant and 533 in sensitive lines, with 186 in common between the two groups (Table S11).Among the latter, many TFs were involved in signal transduction during heat stress (e.g., Heat-Shock Factor 4, HSF4) or in response to hormones, such as abscisic acid, ethylene, and gibberellin.The four TFs predicted only for tolerant lines included two WRKY TFs (WRKY38 and WRKY26); whereas the 347 TFs predicted only for sensitive lines revealed enrichments for biosynthetic and developmental (postembryonic and floral) processes, including ethylene and auxin responses (Figure S6).
We next tested whether any of the identified TFs showed temperature-responsivity in tolerant and sensitive lines.A total of 66 predicted TFs also displayed differential expression in at least one contrast (Table S11, Figure S7).Among those, two families are also involved in heat responses among their numerous functions: the bHLH (seven) and bZIP (eight) TF families.We also found significant TFs involved in responses to auxin, responses to ethylene, multicellular organism development, and flower development.To ease visualisation, expressions of W. ceracea paralogs (orthologs with the same A. thaliana ID) were averaged (89 A. thaliana IDs).Cells are coloured by Z score.Tolerant_Warm = tolerant lines grown at warm temperatures; Sensitive_Warm = sensitive lines grown at warm temperatures; Tolerant_Cool = tolerant lines grown at cool temperatures; Sensitive_Cool = sensitive lines grown at cool temperatures.Labels on the right show respectively W. ceracea gene IDs, A. thaliana gene symbols, and GO Processes.

Discussion
In this study we provided de novo genome and transcriptome assemblies of the Australian endemic alpine herb Wahlenbergia ceracea and functionally annotated thousands of its genes by orthology with Arabidopsis thaliana.We then conducted differential gene expression of thermally tolerant and sensitive lines in response to moderate, prolonged warming to identify candidate genes involved in acclimation and thermal tolerance.Thousands of genes in W. ceracea had significant changes in expression in response to temperature.Overall, we found that tolerant lines downregulated chlorophyll biosynthesis and photosynthesis-related processes, including genes encoding components of the antennae and the electron transport chain when exposed to moderate warming.
Below we discuss the importance of these results for acclimation of photosynthetic thermal tolerance and pinpoint knowledge gaps for further investigation.

Down-regulation of photosynthetic processes in thermally tolerant lines of W. ceracea
Chlorophyll fluorescence is a widely used technique to evaluate the thermostability of the photosystem (PS) II (Geange et al., 2021).Increases in basal chlorophyll fluorescence (F0) observed in gradually heated leaves, corresponding to the Tmax parameter, reflect thylakoid membrane instabilities leading to the discharge of electrons and, therefore, PSII disruption.Less well known are the molecular processes that lead to the acclimation in this parameter observed in many species, including W. ceracea (Notarnicola et al., 2021).
Here, tolerant lines, characterised as having higher Tmax values, responded to warmer conditions by down-regulating genes encoding proteins related to thylakoid membrane organization, pigment metabolism, and light harvesting reactions.These results are confirmed by the analysis of the expression levels of specific photosynthesis-related genes.Specifically, OEC33, ATPC1, and many LHCAs and LHCBs, among other genes encoding components of the PSII or antenna complexes, had lower expressions in tolerant than sensitive lines.A member of the FAD family (FAD6), involved in membrane desaturation in response to temperature, was also downregulated in tolerant lines.These results suggest that part of acclimation to warmer temperatures may include the downregulation of electron transport from PSII.This is considered a photoprotective mechanism to avoid PSI photoinhibition, which impairs photosynthetic assimilation and photoprotection (Zivcak et al., 2015;Jiang et al., 2021).Under heat-stress, the over-reduction of the electron transport chain is compensated by the cyclic electron transport around PSI (Wang et al., 2006;Rumeau et al., 2007;Jiang et al., 2021).Interestingly though, in our study the genes participating in the NAD(P)H dehydrogenase complex assembly, involved in one of two pathways that generate the cyclic electron transport, were down-regulated in tolerant lines but not in sensitive lines.
In tolerant lines we also found a downregulation of the synthesis of chlorophyll and other pigments.This result is partly in agreement with Wang et al. (2018), which conducted wholetranscriptome sequencing on heat-tolerant strains of Pyropia haitanensis exposed to heat stress of several durations, from short-term (a few hours) to long-term (up to six days).However, in Wang et al. (2018) the proportion of downregulated genes related to pigment metabolism decreased with increased heat stress durations.It was proposed that the degradation of photosynthetic pigments and the control of their synthesis may represent physiological mechanisms to reduce the generation of phototoxic intermediates from these compounds during over-excitation (Hörtensteiner, 2006;Wang et al., 2018;Hu et al., 2020).Taken together, these results suggest that photosynthetic thermal tolerance is achieved via conservative reductions of PSII over-excitation and phototoxic damage even at moderately warmer temperatures.
We expected that HSPs would play important roles in the acclimation and thermal tolerance of warmer temperatures but while the 'protein folding' pathway, which includes the activity of HSP chaperones, was significant in overall responses to temperature, we did not observe differences in its regulation between tolerant and sensitive lines.

Prediction of transcription factors reveal possible upstream regulators of thermal acclimation in W. ceracea
Transcription factors (TFs) are proteins that mediate the regulation of gene expression by binding the DNA in regulatory regions upon perception of stress or hormonal signals (Meshi and Iwabuchi, 1995;Haider et al., 2021a).Therefore, predicting the TFs governing expression of DEGs allowed us to evaluate the potential upstream regulators and signals involved in thermal tolerance and acclimation.
We found that many DEGs in both tolerant and sensitive lines are regulated, or participate in, ethylene signalling.Several studies (reviewed in Poór et al., 2021) demonstrated the involvement of ethylene in plant responses to heat stress through regulation of HSFs and HSPs.
For example, in Wu and Yang (2019), rice seedlings treated with a precursor of ethylene had higher thermal tolerance, chlorophyll a content, and higher expression levels of HSF1, HSF2, and antioxidant enzymes under heat stress than seedlings not treated with this precursor.Besides few relevant differences, in our study ethylene responses of tolerant and sensitive lines were mostly regulated by similar members of the ERF family.Future targeted studies of the TFs that do differ between tolerant and sensitive lines may reveal the importance of this hormone in acclimation and thermal tolerance.Generally, the roles of ethylene during regulation of responses to warming warrant further investigations in native species.
It is known that ethylene can induce leaf senescence in plants, both under optimal and high temperature conditions (Li et al., 2013;Kim et al., 2020).Therefore, in annual plants, such as W.
ceracea, and during abiotic stress, this hormone may be involved in systemic responses to increase leaf turn-over and remobilise nutrients from chloroplasts to reproductive organs, seeds, or new leaves (Schippers et al., 2015;Sade et al., 2018).This result agrees with the downregulation of photosynthetic components mentioned above and the higher investment in flowers and growth observed in W. ceracea under warmer temperatures (Notarnicola et al., 2021).
Many DEGs in our study were regulated by other phytohormones, such as cytokinin and jasmonic acid, suggesting these may also be involved in thermal acclimation in W. ceracea.
Cytokinins promote production of some osmolytes (involved in the detoxification of free radicals), stomata opening, and acclimation during heat stress (Dobrá et al., 2015;Escandón et al., 2016;Sharma et al., 2019).During prolonged warming, these hormonal pathways may regulate reactive oxygen species (ROS).Indeed, in our study the enzymatic scavengers commonly found after heatshock, such as APXs and Catalases were not differentially expressed, except for a stromal Ascorbate Peroxidase (sAPX) that was significant in all three contrasts (Table S8 (Zhao et al., 2020;Haider et al., 2021a;Jing et al., 2022) and respond to abscisic acid (Kim et al., 2004;Skubacz et al., 2016;Dixit et al., 2018), a hormone involved in several functions during heat stress, such as induction of antioxidants and regulation of HSPs (Li et al., 2021).Instead, bHLHs were down-regulated in our study, consistently with their functions in cool conditions (Sun et al., 2018).

RNA regulation and post-transcriptional processes in native plants
Post-transcriptional regulation indicates a variety of chemical reactions that control the quantity and quality of mRNA molecules (Mazzucotelli et al., 2008;Guerra et al., 2015).A pivotal role in these processes is played by ribonucleoprotein complexes that are involved in splicing, i.e., the removal of introns and the re-arrangement of exons to generate diverse transcripts from one gene sequence (Mazzucotelli et al., 2008;Guerra et al., 2015), or control mRNA decay to ensure transcript turnover (operated by stress granules and processing bodies; Perea-Resa et al., 2016;Chantarachot and Bailey-Serres, 2018;Kawa et al., 2020).These processes interact with other regulatory mechanisms, such as protein modifications (post-translational regulation), and then feedback back to gene expression adding several layers of regulation to fine-tune responses during thermal acclimation (Miller et al., 2013;Mazzucotelli et al., 2008;Guerra et al., 2015).However, their investigation in native species remains underrepresented.
Evidence from this study revealed that RNA-related processes and ribonucleoprotein complexes may indeed play significant roles in thermal tolerance.Estravis-Barcala et al. (2021) found similar results when comparing Nothofagus pumilio, A. thaliana and Populus tomentosa responses to heat stress.Splicing or RNA and protein modifications induced by heat stress may result in the production of slightly variable protein with altered activity or thermal stability (Reddy, 2007;Staiger and Brown, 2013).Examples of this were found in Rubisco activase (Scafaro et al., 2019;Degen et al., 2021) and HSFA2 (Liu et al., 2013).So far, the involvement of these mechanisms in plant heat tolerance was predominantly studied in Arabidopsis, but our results indicate native species as potential model systems, where the importance of post-translational and post-transcriptional mechanisms for ecological and evolutionary processes can be pinpointed.
Taken together, our results show that responses that occur under a moderate, longer-term heat exposure align with responses to more acute, extreme conditions, although some differences can be pinpointed.However, we note that key genes involved in responses to prolonged warming may go unnoticed if they are not annotated or do not have orthologs with model species.For example, a small fraction (18.6%) of the DEGs in our study did not have orthologs with A. thaliana and may represent novel genes involved in specific responses of W. ceracea or may be pseudogenes or non-coding RNAs (lncRNAs) with post-transcriptional functions.A wider application of these techniques and development of genomic resources for native species will therefore increase our ability to detect differences among species and conditions.Ultimately, studying new genes and/or gene variants in more non-model, natural species will enable us to predict responses to climate change and to engineer more resistant crop varieties.

Figure 2 :
Figure 2: Tolerant lines exhibited stronger enrichment of photosynthesis-related terms among down-regulated DEGs; sensitive lines had stronger enrichment of post-transcriptional processes among up-regulated DEGs.Bar charts of the significant enriched terms: (a) tolerant lines, upregulated terms; (b) sensitive lines, upregulated terms; (c) tolerant lines, downregulated terms; and (d) sensitive lines, downregulated terms.Significant terms were clustered by Metascape; the terms with the highest P-values in each cluster were chosen as representatives.Significant terms are from GO Processes ('GO:…') and KEGG Pathways ('ath…') catalogues.Bars are filled by log(Q-value) calculated for the individual DEG lists (tolerant or sensitive), where Q-value is the adjusted P-value; more negative values (darker bars) correspond to smaller Q-values.Terms are ordered based on their overall Q-values in the combined DEG lists (tolerant + sensitive) in descending order with more significant terms on top.White bars = non-significant terms.Xaxis = fold enrichment.

Figure 3 :
Figure3: Heat map of the average expression of the photosynthesis-associated DEGs shows that most of them were more down-regulated in tolerant lines under warming.To ease visualisation, expressions of W. ceracea paralogs (orthologs with the same A. thaliana ID) were averaged (89 A. thaliana IDs).Cells are coloured by Z score.Tolerant_Warm = tolerant lines grown at warm temperatures; Sensitive_Warm = sensitive lines grown at warm temperatures; Tolerant_Cool = tolerant lines grown at cool temperatures; Sensitive_Cool = sensitive lines grown at cool temperatures.Labels on the right show respectively W. ceracea gene IDs, A. thaliana gene symbols, and GO Processes.

Table 1 :
Differential gene expression contrasts performed and their biological interpretation.DEGs = differentially expressed genes; Tolerant_warm = tolerant lines grown at warm temperatures; Sensitive_warm = sensitive lines grown at warm temperatures; Tolerant_cool = tolerant lines grown at cool temperatures; Sensitive_cool = sensitive lines grown at cool temperatures.