Deciphering the transcriptomic regulation of heat stress responses in Nothofagus pumilio

Global warming is predicted to exert negative impacts on plant growth due to the damaging effect of high temperatures on plant physiology. Revealing the genetic architecture underlying the heat stress response is therefore crucial for the development of conservation strategies, and for breeding heat-resistant plant genotypes. Here we investigated the transcriptional changes induced by heat in Nothofagus pumilio, an emblematic tree species of the sub-Antarctic forests of South America. Through the performance of RNA-seq of leaves of plants exposed to 20°C (control) or 34°C (heat shock), we generated the first transcriptomic resource for the species. We also studied the changes in protein-coding transcripts expression in response to heat. We found 5,214 contigs differentially expressed between temperatures. The heat treatment resulted in a down-regulation of genes related to photosynthesis and carbon metabolism, whereas secondary metabolism, protein re-folding and response to stress were up-regulated. Moreover, several transcription factor families like WRKY or ERF were promoted by heat, alongside spliceosome machinery and hormone signaling pathways. Through a comparative analysis of gene regulation in response to heat in Arabidopsis thaliana, Populus tomentosa and N. pumilio we provide evidence of the existence of shared molecular features of heat stress responses across angiosperms, and identify genes of potential biotechnological application.


Introduction
show a time of day-dependent response, and it was demonstrated that both diurnal and 86 circadian regulation of the transcriptome impact experimental interpretation of the heat 87 stress response [27]. With the aim to detect genes regulated by heat stress in N. 88 pumilio, and reduce the aforementioned possible diurnal (photocycles-driven) and time 89 of the day-dependent bias in the interpretation of the heat stress experiments, we used 90 the following protocol. Plants were grown for 10 days in growth chambers (SCE Sequencing was performed using the ionTorrent Proton System (Life Technologies), 122 in a total of three runs (with three, three, and two libraries, respectively) in order to 123 ensure approximately 25 million reads per library, which was shown to be sufficient to 124 detect more than 90% of genes in eukaryotes [29]. 125 Datasets processing and assembly 126 Reads were quality-checked with FastQC [30] and trimmed with Trimmomatic [31]   Trimmed reads were assembled using SPAdes [33] (version 3.11.0; parameters: --rna 131 --iontorrent -k67 --ss-fr). The final k -mer value of 67 was chosen after several assemblies 132 with different k -mer values (five in total, between 21 and 77). The assembled contigs 133 that overlapped considerably were expanded, and highly redundant contigs were 134 eliminated. 135 After assembly, trimmed reads were mapped back to the assembly using STAR [34]  contigs were aligned to the database (A. thaliana proteome) using BLAT [37], and a file 156 with a single best-hit annotation for each successful contig was generated. The  Reads from all libraries were quantified against the reference assembly using Salmon, 162 version 0.8.1 [39]. After quantification, a tab-delimited file containing the unnormalized 163 expression level for each contig in each of the eight libraries was put together. For 164 differential expression analysis, DESeq2 [40] was used in an R environment, with default 165 models and parameters. The two temperatures (20 • C and 34 • C) were contrasted, 166 taking the different moments of the day as biological replicates; that is, four biological 167 replicates for each temperature were compared. This protocol allowed us to reduce the 168 bias of the time of the day on the interpretation of our experiments, focusing our study 169 in the detection of genes that were particularly induced by heat stress. Contigs with an 170 FDR<0.05 were considered as differentially expressed between temperatures.
Functional enrichment analysis assembly, we were able to perform GO terms and metabolic pathways enrichment in 174 contigs over-and under-expressed in response to temperature. For GO term enrichment 175 analysis, PANTHER version 14.1 [41] was used via its implementation in the TAIR 176 (The Arabidopsis Information Resource; https://www.arabidopsis.org/) website.

177
For KEGG (Kyoto Encyclopedia of Genes and Genomes; [42]) metabolic pathways 178 enrichment analysis, KOBAS 3.0 online tool was used [43]. In both cases, the 179 background dataset were all A. thaliana identifiers present in our assembly's annotation. 180 Visualization and clustering of over-represented GO terms was performed with 181 REVIGO [44].  In order to inspect shared molecular components that are up-or down-regulated in 191 response to heat stress among plant species, we used expression data from A. 192 thaliana [47] and Populus tomentosa [18]. These papers were selected among those 193 published in recent years because they feature a complete, publicly available set of gene 194 annotations and differential expression statistical results. The P. tomentosa study 195 consisted in RNA-seq experiments where contigs were annotated against P. trichocarpa, 196 a related species with a sequenced genome, which in turn was annotated against A.

197
thaliana. Thus, we were able to obtain corresponding A. thaliana IDs for annotated 198 contigs for all three species that could be intersected and provided us with a list of 199 shared genes in the three species. Among these, each species had a set of differentially 200 expressed genes, which were also intersected and subjected to GO term enrichment and 201 transcription factor regulation prediction as described above for N. pumilio.
Primer design and quantitative RT-PCR validation 203 In order to evaluate the accuracy of our transcriptome data, a total of 13 (eight 204 up-regulated and five repressed in response to high temperature) genes were selected to 205 carry out a qRT-PCR analysis. We chose a group of contigs that allowed to test a wide 206 range of expression (from intermediate to high expression; between 5 and 45 Transcripts 207 Per Million averaged across conditions) and fold change (Log2 fold change between -11 208 and 8) in our RNA-seq data. Primers were designed with Primer-BLAST [48] (S1 209   Table).   Table 1). The read 226 throughput and average length were in accordance with the device specifications [52].

227
Raw reads were trimmed to eliminate low-quality ends, and low-quality sequences 228 (Q<20) were removed. This procedure allowed us to increase the overall read quality at 229 the expense of shorter overall read length (  Table).  Table). This reference transcriptome was then used for differential expression and and accounted for 1,265 and 883 unique protein IDs, respectively (S13 Table and S14   260   Table).  Table) were directly related to photosynthesis, for 269 example "Photosynthesis", "Photosynthesis -antenna proteins", "Carotenoid in genes promoted by 34 • C ( Fig 1B) were specific to the biosynthesis of stress-related metabolites families such as flavonoids, phenylpropanoids, mono-, sesqui-and 280 tri-terpenoids, and some groups of alkaloids. Pathways related to translation and 281 protein processing were also triggered at 34 • C, as indicated by several enriched 282 pathways such as "Ribosome", "Protein processing in endoplasmic reticulum", 283 "Spliceosome" or "RNA transport" (Fig 1B and C, S5 Table). Overall, a total of 78 GO terms were enriched in genes repressed at 34 • C (S6 Table) 285 and 71 GO terms were enriched in genes induced by 34 • C (S7 Table). Semantic 286 reduction and clustering of enriched GO biological processes show that 287 "photosynthesis", "glucose metabolism" and "generation of precursor metabolites and 288 energy" were the main processes repressed by heat, whereas at 34 • C the response to 289 various stress signals such as "response to chemical" or "secondary metabolism" and 290 "protein folding / refolding" were highly significant (Fig 2). Moreover, the most enriched 291 GO terms in genes repressed at 34 • C in all three branches of the ontology (Biological  Table and S7 Table). Tables 2 and 3    Hormones play a fundamental role in plant stress responses, and the function of 318 particular hormones and their crosstalk differ among tree species [3]. The homolog of 319 AHP5 (Histidine-containing phosphotransfer protein 5), an important two-component 320 mediator between cytokinin sensing and its response regulators, was repressed at 34 • C 321 (S13 Table). In addition, ARFs TF family was down-regulated at 34 • C (S8 Table). (S14 Table). These ABA phosphatases, which show high homology to the ABA 327 phosphatases At4g26080, At5g59220, At1g07430 and At3g11410 of A. thaliana, are part 328 of the KEGG pathway "MAPK signaling pathway -plant", over-represented in genes 329 induced by 34 • C (Fig 1B).

330
Comparative heat stress responses between species

331
The comparative analysis of heat stress response of N. pumilio, A. thaliana and P. 332 tomentosa) resulted in 68 genes that were significantly more expressed at high 333 temperature in all three species (S10 Table), many of them belonging to the 334 aforementioned up-regulated groups in N. pumilio, such as HSPs, LEAs (Table 2) the DnaJ protein ERDJ3A or the DNAJ protein P58IPK homolog, that contribute to 342 the protection of cells to endoplasmic reticulum stress [54].

343
GO enrichment analysis showed that the main processes shared among species in 344 response to heat were those related to protein misfolding and refolding, apart from 345 general and specific stress terms (S11 Table and  high temperature (S12 Table). Of those, 29 (55%) were ERFs, demonstrating the 348 relevant role of ethylene in the response to high temperature stress in these species.

349
Validation of RNA-seq data with quantitative RT-PCR

350
To verify the validity of our RNA-seq differential expression results, we analysed the These results show the high reliability of the RNA-seq data. ribosome, peptide biosynthesis, and translation (Fig 1, Fig 2). A well-known effect of 390 abiotic stress is the production of Reactive Oxygen Species (ROS), which can oxidize 391 biomolecules and set off cell death [60]. Many plant secondary metabolites have 392 antioxidant properties, and their production is significantly increased by abiotic 393 stress [60,61]. In our study, genes involved in the biosynthesis of many antioxidant 394 metabolites families were found to be triggered by high temperature, namely phenylpropanoids, flavonoids, mono-, sesqui-and tri-terpenoids, and tropane, piperidine 396 and pyridine alkaloids. Moreover, it has been shown in trees that MAPK cascades 397 promote antioxidant responses [62], and the plant MAPK signaling pathway was 398 enriched in genes more expressed at 34 • C in N. pumilio (Fig 1B and S5 Table). 399 Our analysis indicated that the response to misfolded or topologically incorrect 400 proteins was up-regulated by heat stress. In plants, HSPs and other chaperones bind to 401 misfolded proteins, which are in turn ubiquitinated and directed to the proteasome for 402 their degradation [63,64]. Several chaperones (including HSPs) and ubiquitin-ligases 403 were found to be significantly more expressed at 34 • C than at 20 • C (Tables 2 and 3), 404 indicating the importance of these processes in the response to high temperature stress. 405 In concordance, forestry species such as Quercus lobata, Pseudotsuga menziesii and

406
Prunus persica show over-expression of chaperones and ubiquitin-ligase proteins in 407 response to different abiotic stresses [13,56,65], indicating that the induction of protein 408 re-folding and ubiquitination followed by degradation by the proteasome pathway 409 constitutes a relevant molecular strategy that allows trees to cope with adverse abiotic 410 conditions.

411
Transcription Factors (TFs) are known to play important roles in the transcriptional 412 regulation of stress responses, and their involvement in many biotic and abiotic stresses 413 across plant species has been extensively reviewed [53,66,67]. In N. pumilio, TFs 414 belonging to families such as WRKY, WOX, LBD and NAC were positively regulated in 415 response to high temperature (S8 Table). In concordance, previous reports show that 416 NAC TFs play roles in numerous biotic and abiotic stresses including heat [68,69], and 417 whereas WOX TFs promotes the response to abiotic stresses such as drought and 418 salinity in Brassica napus [70,71], LBD TFs were suggested to play roles in the response 419 to cold in Broussonetia papyrifera [72]. The WRKY gene family is one of the largest TF 420 family in plants, playing roles in the regulation of a broad range of physiological and 421 developmental processes [73], including the response to biotic and abiotic stress [74,75]. 422 It is interesting to note that most of the WRKY TFs identified in N. pumilio involved in the response to salt stress in poplar trees [80]. Moreover, WRKY18, 48 and 427 over-representation of KEGG pathways and GO terms related to response to oxidative 429 stress at 34 • C (Fig 1, Fig 2), raising the hypothesis that ROS may induce the 430 expression of a sub-set of WRKY TFs during the response to high temperature in N.  Arabidopsis [83], was strongly repressed by heat stress in the N. pumilio transcriptome 436 (S13 Table).

437
In relation to hormonal signaling, ethylene is an important plant hormone which is 438 known to be involved in stress responses [84]. Several members of the ERF family were 439 over-expressed under heat stress in N. pumilio (S8 Table and S14 Table), and genes 440 up-regulated at 34 • C showed an enrichment of ERF targets (S9 Table). Moreover, the 441 homolog of EIN3, a master regulator of ethylene signaling [85], together with several of 442 its targets [86], were over-expressed in heat-treated plants (S14 Table), indicating that 443 the ethylene pathway is activated at 34 • C. Additionally, our data suggest that ABA 444 signaling and response constituted another hormonal pathway up-regulated by heat.

445
This is supported by the over-expression of several ABA-responsive genes such as those 446 described in [53], including dehydrins, LEAs and protein phosphatases of the clade A in 447 plants exposed to high temperature (Table 2 and S14 Table). Furthermore, genes 448 over-expressed at 34 • C showed enriched targets of TFs related to the regulation of ABA 449 signaling (S9 Table). In contrast, our data indicates that auxin signaling and 450 re-localization is repressed in response to heat in N. pumilio. This is evidenced by the 451 fact that ARFs, which are relevant components of the auxin signaling pathway [87], 452 were repressed by high temperature (S8 Table). In addition, auxin-efflux ABC 453 transporters were down-regulated under heat stress (Fig 1A). Finally, RVE1 (Reveille 1), 454 a MYB-like TF which links the circadian clock with the auxin signaling pathway [88], combined analysis of genes over-expressed under heat stress in N. pumilio, A. thaliana 460 and P. tomentosa allowed us to identify a core of shared responses to high temperature, 461 mostly related to protein misfolding and chaperone activity, with the over-expression of 462 more than ten HSPs, one LEA and two DnaJ protein genes (S10 Table,  to AS [90]. In our study, the "Spliceosome" KEGG Pathway was significantly enriched 466 among genes up-regulated at 34 • C in N. pumilio (Fig 1B and S5 Table), and several 467 genes involved in pre-mRNA splicing were shared between N. pumilio, A. thaliana and 468 P. tomentosa at high temperature, including several DEAD-box ATP-dependent RNA 469 helicases (S10 Table). These results support the fact that AS constitutes an important 470 mechanism in plant response to abiotic stress and highlight the potential relevance of a 471 subset of genes associated with the splicing machinery in the response to heat stress 472 across angiosperms.

473
Regarding hormone signaling, many ERFs were found to have enriched targets 474 among genes over-expressed at high temperature in the three species (S12 Table), 475 indicating the relevance of ethylene and ERF TFs in the response to heat, and further 476 supporting the reported results in N. pumilio. Apart from ERFs, our analysis allowed 477 us to identify common targets or relevant TF families already discussed such as WRKY 478 and NAC (S12 Table), and one of the TFs with most significantly enriched targets, and 479 the single most enriched considering only N. pumilio (S9 Table)  be differentially expressed in response to cold, heat and drought treatments in rice [91]. 482 Finally, the zinc finger protein ZAT10, which constitutes a transcriptional repressor 483 involved in abiotic stress responses [92], was over-expressed under heat stress in the 484 three species (S10 Table), and A. thaliana, P. tomentosa and N. pumilio transcriptomes 485 of heat-treated plants show an enrichment of ZAT10 targets (S12 Table). This suggests 486 that the ZAT10 regulon constitutes a relevant regulatory module during heat stress 487 responses in angiosperms. All these results suggest a strong shared core of