Gene expression analysis of Cyanophora paradoxa reveals conserved abiotic stress responses between basal algae and flowering plants

The glaucophyte Cyanophora paradoxa represents the most basal member of the Archaeplastida kingdom, however the function and expression of most of its genes are unknown. This information is needed to uncover how functional gene modules, i.e. groups of genes performing a given function, evolved in the plant kingdom. We have generated a gene expression atlas capturing responses of Cyanophora to various abiotic stresses. This data was included in the CoNekT-Plants database, enabling comparative transcriptomic analyses across two algae and six land plants. We demonstrate how the database can be used to study gene expression, co-expression networks and gene function in Cyanophora, and how conserved transcriptional programs can be identified. We identified gene modules involved in phycobilisome biosynthesis, response to high light and cell division. While we observed no correlation between the number of differentially expressed genes and the impact on growth of Cyanophora, we found that the response to stress involves a conserved, kingdom-wide transcriptional reprogramming, which is activated upon most stresses in algae and land plants. The Cyanophora stress gene expression atlas and the tools found in https://conekt.plant.tools/ database provide a useful resource to reveal functionally related genes and stress responses in the plant kingdom.

• The glaucophyte Cyanophora paradoxa represents the most basal member of the 23 Archaeplastida kingdom, however the function and expression of most of its genes are 24 unknown. This information is needed to uncover how functional gene modules, i.e. groups of 25 genes performing a given function, evolved in the plant kingdom. 26 • We have generated a gene expression atlas capturing responses of Cyanophora to 27 various abiotic stresses. This data was included in the CoNekT-Plants database, enabling 28 comparative transcriptomic analyses across two algae and six land plants. 29 • We demonstrate how the database can be used to study gene expression, co-expression 30 networks and gene function in Cyanophora, and how conserved transcriptional programs can 31 be identified. We identified gene modules involved in phycobilisome biosynthesis, response to 32 high light and cell division. While we observed no correlation between the number of 33 differentially expressed genes and the impact on growth of Cyanophora, we found that the 34 response to stress involves a conserved, kingdom-wide transcriptional reprogramming, which 35 is activated upon most stresses in algae and land plants. 36 tend to have similar expression patterns across organs, developmental stages and biotic as well 72 as abiotic perturbations (Usadel et al., 2009). Co-expression analyses are widely used 73 (Radivojac et al., 2013;Rhee & Mutwil, 2014;Jiang et al., 2016), and have been applied to 74 successfully identify genes involved in plant viability (Mutwil et al., 2010), seed germination 75 (Bassel et al., 2011), shade avoidance (Jiménez-Gómez et al., 2010), cyclic electron flow 76 visualizing the starch levels, the cells were stained with 20% Lugol solution and analysed by 133 microscopy (Fig. 1c). For the optical density (OD) measurements, 20 mL of culture were μ L iodine reagent (1% KI + 0.1% I 2 ) was 139 added right before measuring the absorbance at 680 nm (Avidan et al., 2015). Three 140 independent cultures with different cell densities were used for this analysis: i) ~0.5 * 10 6 141 cells/mL, ii) ~10 6 cells/mL and iii) ~2*10 6 cells/mL. Three technical replicates were 142 performed. 143 6 144

RNA preparation and sequencing 145
The RNA was extracted using Spectrum™ Plant Total RNA Kit (Sigma-Aldrich) according to 146 the manufacturer's instructions. The integrity of the RNA was assessed using RNA nano chip 147 on Agilent Bioanalyzer 2100. The libraries were prepared from total RNA using polyA 148 enrichment and sequenced using Illumina-HiSeq2500/4000 at Beijing Genomics Institute and 149 Max Planck-Genome-centre in Cologne. Three RNA isolations were done for each sample. 150 151

Analysis of RNAseq data and microarrays data 152
The reads were trimmed, mapped, counted, and TPM-normalized using the LSTrAP pipeline 153 (Proost et al., 2017). The genome used for mapping the reads is C. paradoxa v1.0 (Price et al., 154 2019). More than 81% of the reads mapped to the genome, and on average 88% of the reads 155 mapped to the coding sequences (Table S2). The principal component analysis revealed that 156 one of the three replicates of salt stress did not cluster with the other two replicates. For this 157 reason, we excluded the outlier sample and sequenced two additional biological replicates. 158 The publicly available data was processed as follows: raw expression data for 159 Arabidopsis thaliana, Chlamydomonas reinhardtii, and Synechocystis sp. PCC 6803 (Table  160 S3) were downloaded from ArrayExpress. The raw fastq files were processed using the 161 LSTrAP pipeline and the reads were mapped to A. thaliana TAIR10 and C. reinhardtii v5.5. 162 Raw counts from the RNA-seq experiments were used with the R package DESeq2 (Love et 163 al., 2014) to identify differentially expressed genes (DEGs). The various types of microarrays 164 were processed based on the manufacturer. Affymetrix cel files were processed with the R 165 package affy (Gautier et al., 2004). Agilent raw files were processed with the R package 166 limma (Ritchie et al., 2015). Differentially expressed genes for all microarray experiments 167 were identified by using the R package limma. For further analyses, we only considered the 168 genes which showed an adjusted p-value < 0.05 and a -1>log 2 fold>1 as DEGs (Table S4-S7). The normalized TPM matrix containing expression values from diurnal cycle (Ferrari et al.,173 2019) and the expression matrix of the newly generated stress samples (Table S8)  The Heuristic Cluster Chiseling Algorithm (HCCA) clusters of C. paradoxa were downloaded 184 from CoNekT-Plants (Table S9). For each cluster we counted the number of times a specific 185 biological process was represented. This observed distribution of biological processes was 186 compared to a permuted distribution obtained by sampling an equal number of genes from the 187 total pool of C. paradoxa genes for 10,000 permutations. The empirical p-values obtained 188 were false-discovery rate (FDR) corrected by the Benjamini-Hochberg procedure (Benjamini 189 & Hochberg, 1995). 190

Gene function analysis of abiotic stress 192
To investigate how the stress conditions perturbed certain biological processes (MapMan 193 bins), we first analyzed the fraction of differentially expressed genes (DEGs) assigned to each 194 bin. We assigned the fractions into strongly affected (≥50% of genes in the bin are DEGs), 195 mildly affected (20% ≤ DEGs < 50%) or weakly affected (DEGs < 20%). To visualize if the 196 DEGs found in the bin are preferentially up-or down-regulated we used the formula: 197 The formula returns 1 when all genes in the bin are upregulated, 0 when there is an equal 198 number of up-and down-regulated genes and -1 when all genes in the bin are downregulated. 199

Inferring similarities of stress responses 201
We used the Jaccard Index (JI) to assess the similarity of stress responses between two 202 experiments. To this end, we first identified the gene families corresponding to the DEGs, 203 which were used to calculate the Jaccard index for two compared stresses, X and Y. The 204 observed JI between stresses X and Y is: 205 To test whether a given JI obs(X,Y) is larger than expected by chance, i.e. two stresses show 207 similar transcriptomic response, we performed a permutation test where we drew a number of 208 genes equal to the amount of DEGs present in a specific experiment (DEGs_r) and calculated 209 the expected Jaccard indexes for each experiment pair (JI perm(X,Y) ). The procedure was 210 performed 10,000 times, and the empirical p-values obtained, calculated independently for up-211 and downregulated families, were FDR corrected. The OrthoFinder gene families used in this 212 analysis were obtained from (Ferrari et al., 2019). 213 214

Identification of stress-responsive gene families 215
We counted the number of experiments when a given family contained DEGs, and calculated 216 the empirical p-value by estimating if this number was larger than the corresponding counts 217 obtained from shuffled gene-DEG assignments. The shuffling was done 10,000 times, and the 218 empirical p-values were FDR corrected. 219 220

Phylostratigraphic enrichment analysis of DEGs 221
To test whether there is a correlation between the age (phylostrata) (Domazet-Loso et al., 222 2007) of a gene family and the transcriptomic responses to stress, we calculated the fraction of 223 DEGs found in each phylostratum. Every fraction was then tested for significance by 224 estimating whether the observed fraction was larger than fractions of shuffled gene-phylostrata 225 assignments. The shuffling was done 10,000 times, and the empirical p-values were FDR 226 corrected. The gene-phylostrata assignment was retrieved from (Ferrari et al., 2019). 227 228

Data availability 229
The raw sequencing data is available from EBI accession number E-MTAB-7822. 230 231 Results 232

Generation of C. paradoxa expression atlas 233
To be able to assign gene functions via co-expression analyses we induced broad changes in 234 gene expression in C. paradoxa by subjecting the liquid culture to four nutritional, one 235 osmotic, one salt and four environmental stresses. The intensities of the different stresses were 236 modulated to reduce the growth rate without stopping it entirely. To this end, culture growth 237 was monitored by cell counting from the time when the stress was first applied (time 0 hours) 238 up to several days after the harvest time for gene expression analyses (Fig. 1a, b, and d). The 239 nutritional stress category involves the depletion of macro or micronutrient to the growth 240 medium (nitrogen, sulphate, phosphate, micronutrient solution), while the salt and osmotic 241 stresses were induced by addition of 75 mM NaCl or 100 mM mannitol, respectively. Harvest 242 time for all nutritional, salt and osmotic stress conditions was on day 4 when a pronounced 243 reduction of growth was observed for most stresses compared to the control (Fig. 1a). The 244 presence of NaCl in the medium had the strongest negative effect on growth, followed by 245 presence of mannitol, depletion of phosphate and depletion of TMS. The effect of sulphate 246 depletion was only visible on day 5 (Fig. 1a). For the environmental stresses we subjected the 247 cultures to heat (37°C), cold (4°C), high light (150 μ E) with a harvest time at 12 h for each 248 stress condition ( Fig. 1b and 1d). Cultures grown for 72 hours in darkness showed an expected 249 absence of starch granules as compared to the control (Fig. 1c, Fig. S1). The environmental 250 stresses resulted in less pronounced negative growth effects (Fig. 1b, d), however on day 2.5 251 and day 3.5 after the environmental stress was applied, the effects were more clearly visible at 252 the level of growth rate (Fig. 1b, d). 253 To evaluate how gene expression in Cyanophora was affected by the different stresses, 254 RNA was isolated from the harvested cultures and sequenced in triplicates to produce at least 255 15 million reads per sample (Table S2) To identify the significantly differentially expressed genes (DEGs, adjusted p-value < 261 0.05, -1>log 2 fold>1), we compared the raw counts of the stress conditions to the respective 262 control. The impact of the stress on the transcriptome was highly variable, ranging from 270 263 DEGs detected after 12 hours of high light treatment to 7628 DEGs after 72 hours of darkness 264 (Table S5). Surprisingly, we did not observe any significant correlation between the number of 265 DEGs and the magnitude of effect the stress had on the relative growth of the culture (Fig. 1e). 266 Cold and sulphate deprivation (gray triangle and red dot, respectively) showed mild inhibition 267 of growth, but a much higher proportion of DEGs compared to treatments with a severe 268 growth phenotype, such as salt and mannitol treatments (Fig. 1e, turquoise and orange dots 269 respectively). Consequently, we observed no significant correlation between DEGs and 270 relative growth rate in Cyanophora (r = 0.278, p-value = 0.437). 271 272

Integration of Cyanophora transcriptome into CoNekT database 273
Biological networks are characterized by their scale-free topology nature, which results in few 274 genes being connected (correlated) to many genes, while the majority of genes have only a few 275 connections (Barabási & Bonabeau, 2003). This architecture is presumed to ensure stability in 276 the case of perturbation, as the network topology remains unaffected when certain genes 277 mutate ( Barabási & Oltvai, 2004). To test if the produced C. paradoxa data also follows an 278 expected scale-free network, we calculated the Pearson Correlation Coefficient (PCC) for 279 every gene pair, with a threshold of 0.9 (Usadel et al., 2009), and counted the number of times 280 a given gene is co-expressed with other genes at this threshold (node degree). The points of the 281 resulting power law plot formed a line with a negative slope, showing a negative correlation 282 between node frequency (i.e., number of genes with a certain number of connections) and 283 node degree (i.e., number of connections per gene). Thus, this result confirms the scale-free 284 topology of the network and indicates that our expression data produces biologically relevant 285 topology (Fig. 2a). 286 In order to easily browse the expression profiles of genes under different stress 287 conditions, to analyze the Cyanophora co-expression network (Table S10) and functional gene 288 modules, and to compare the gene modules to higher species of the Archaeplastida kingdom, 289 we added our gene expression data, together with publicly available diurnal samples for 290 Cyanophora (Ferrari et al., 2019) to the CoNekT-Plants platform (Proost & Mutwil, 2018). 291 CoNekT-Plants is a user-friendly web-tool containing functional and expression data for 8 292 species including a glaucophyte (the newly added Cyanophora paradoxa), chlorophyte 293 (Chlamydomonas reinhardtii), a gymnosperm (Picea abies), two monocots (Oryza sativa, Zea 294 mays) and three dicots (Vitis vinifera, Arabidopsis thaliana and Solanum lycopersicum), and 295 allows comparative genomic and transcriptomic network analysis for these species. CoNekT-296 plants allows the users to view expression profiles and co-expression networks of their genes 297 of interest, but also performs more sophisticated analyses, such as identification of functional 298 gene modules involved in a biological process of interest, comparison of gene expression 299 across different species, identification of genes highly expressed in a given stress or organ, and 300 others (Proost & Mutwil, 2018). 301 To exemplify the applicability of the tool for Cyanophora, we explored the expression 302 profile and co-expression network of CpcK2 (this was done by entering 303 and is existing to optimize the expression and activity of genes during the day/night cycle 317 (Bell-Pedersen et al., 2005). 318 We explored the co-expression network associated with the CpcK2 (by clicking on the 319 graph icon under the 'Co-expression networks/Neighborhood' panel), which showed genes 320 directly co-expressed (direct neighbors) with the query gene. In the network we identified 321 additional linker protein encoded genes (CpcK1, CpcG1, CpcG2, CpcD, ApcC1, ApcC2) 322 Thus, the database allows to rapidly identify functionally related genes in Cyanophora. 328 329

Functional analysis of network clusters 330
Functionally related genes form highly connected clusters in the co-expression network, and 331 these clusters can be identified and studied to unravel the functional gene modules of an 332 organism (Mutwil et al., 2010;Rhee & Mutwil, 2014;Aoki et al., 2015) To identify the 333 functional gene modules of Cyanophora, we used the Heuristic Cluster Chiseling Algorithm 334 (HCCA) (Mutwil et al., 2010), and obtained 319 clusters of co-expressed genes (Table S9). To 335 reveal their function, we performed an enrichment analysis and identified 114 clusters 336 significantly enriched for at least one biological process (Fig. 3a, Fig. S3 Nucleotide metabolism were enriched in multiple clusters, suggesting a more complex 342 regulation of these processes involving several gene modules (Fig. S3). 343 CoNekT-Plants shows the average expression of all genes belonging to a given cluster. 344 As an example, we explored cluster 54 (https://conekt.sbs.ntu.edu.sg/cluster/view/45), which 12 was significantly enriched for the process "Photosynthesis" (Fig. 3a, indicated by a gray 346 rectangle), and observed that the cluster expression is the highest for the sample 347 "Environmental stress -High light -150 μ E". In contrast, "Environmental stress -Prolonged 348 darkness -72h" showed the lowest expression in this cluster, supporting the evidence of 349 association with photosynthesis-related processes (Fig. 3b). This observation suggests that 350 cluster 54 contains genes that are highly responsive to the amount of light perceived by 351

Cyanophora. 352
We further explored the network associated with cluster 54 and by looking at the labels 353 obtained from InterPro annotation, we identified five major classes of genes: oxidoreductases, 354 transporters, photosynthesis machinery related proteins, peptidases, and DNA/RNA repair and 355 modification proteins (Fig. 3c). Additionally, Gene Ontology enrichment found on the cluster 356 page for cluster 54 (GO enrichment is found under the expression profile), showed a 357 significant enrichment of terms related to oxidoreduction ("oxidation-reduction process" 358 GO:0055114, "oxidoreductase activity" GO:0016491, "oxidoreductase activity, acting on 359 single donors with incorporation of molecular oxygen" GO:0016701). These results suggest 360 that this light responsive cluster is likely to be active during high light stress and is involved in 361 the cellular repair of biological components damaged by high light. Thus, this example 362 underlines the applicability of co-expression analyses on uncovering functional gene modules. 363 364

Comparative co-expression analysis reveals conserved cell division program 365
For another case study illustrating comparative co-expression analyses in Cyanophora we 366 selected the cell cycle and its regulation, which have been extensively studied in many model 367 organisms. Notably, the many components involved in DNA replication and chromosome 368 separation are highly conserved across clades and kingdoms of life (Stuart et al., 2003;Yu et 369 al., 2003). However, given the complexity of multicellular organisms, certain aspects of cell (https://conekt.sbs.ntu.edu.sg/tree/view/5112) revealed that the gene has orthologous genes in 384 C. reinhardtii, S. moellendorfii P. abies, O. sativa, Z. mays, V. vinifera, A. thaliana, and S. 385 lycopersicum, and that these orthologs are expressed in various tissues in land plants (Fig. 4a). 386 CoNekT uses the Expression Context Conservation (ECC) and label co-occurrence 387 values to identify conserved co-expressed gene modules (Movahedi et al., 2011;Proost & 388 Mutwil, 2018

Comparative analysis of stress responses in algae and flowering plants 413
Comparative transcriptomic studies across species have helped us considerably to understand 414 which biological pathways are conserved and how they have evolved over time (Mutwil et al., 415 2011;Movahedi et al., 2012;Hansen et al., 2014;Ruprecht et al., 2017a). For example, we 416 could show that diurnal transcriptional programs are conserved over a huge time span of 1 417 billion years (Ferrari et al., 2019). The co-expression networks are also conserved across 418 several species (Ruprecht et al., 2011(Ruprecht et al., , 2016(Ruprecht et al., , 2017aMutwil et al., 2011). In the present study 419 we aimed to investigate if the same is true for transcriptomic responses to a wide range of 420 abiotic stresses. To this aim, we retrieved publicly available gene expression data capturing 421 similar stresses as studied in C. paradoxa (Table S3)  is available for C. paradoxa and A. thaliana. We also included our data on mannitol stress for 426

C. paradoxa. 427
We first investigated how many genes change their expression profile upon exposure 428 to different stresses. The fraction of stress-specific DEGs is overall below 20% across species, 429 except for prolonged darkness, which results in about 26% and 31% of DEGs in Arabidopsis 430 and Cyanophora, respectively (Fig. 5a, gray bars, Table S4-S7). Other stresses that caused 431 strong transcriptional responses were nitrogen deprivation and temperature stresses. 432 Interestingly, Synechocystis shows an overall large transcriptomic response with many DEGs 433 found for most stresses, suggesting that cyanobacteria might respond more strongly to stresses. 434 We further observed that in most cases the proportion of upregulated and downregulated genes 435 was equal (Fig. 5b,  We then explored how the different biological processes were affected by these 440 stresses. To test that, we analyzed the fraction of DEGs in each of the MapMan ontology bins 441 (Fig. 5c). Prolonged darkness had the strongest impact on the transcriptome, especially for C. 442 paradoxa, where more than one third of all biological processes (11 out of 29) showed more 443 than 50% of DEGs (Fig. 5c). The majority of these processes shows a mild prevalence of 444 downregulated genes, suggesting a response to carbon starvation. Interestingly, we observed a 445 general tendency of downregulation of genes involved in photosynthesis for all four species 446 and across all stress conditions. 447 To test whether stress responses are conserved across species, we measured the 448 similarity of responses of gene families in the different experiments. As a measure of 449 similarity we determined whether a pair of stresses contains more common up-or down-450 regulated gene families than expected by chance. We found that the down-regulative responses 451 (blue cells) are more frequent than the up-regulative responses (red cells, Fig. 6a). The high 452 light treatment shows very little significant similarity to the other stresses, especially in C. 453 reinhardtii and C. paradoxa. Conversely, the other conditions showed a significantly similar 454 down-regulatory responses across the stress types and species. For example, the nitrogen 455 depletion response is similar across the four analyzed species, but also similar to the other 456 nutrient depletion conditions, heat, and cold responses (Fig. 6a,

blue cells). Overall we 457
conclude that the down-regulatory responses are highly conserved, regardless of the type of 458 stress and species. 459 Next, we set to investigate which gene families tend to significantly respond to the 460 stresses ( Figure 6b). We observed that gene families involved in photosynthesis, tetrapyrrole 461 biosynthesis, and glycolysis, showed a predominant downregulation response. Other 462 processes, such as nutrient uptake, cellular respiration, amino acid, and protein biosynthesis To investigate whether the age of a gene family is correlated to its stress response, we 476 analyzed the fraction of genes that is either up-or downregulated in each phylostratum (i.e., 477 the relative age of the gene family). We arranged the phylostrata from oldest (prokaryotic) to 478 youngest (species-specific gene families) and observed that the response to stress is evenly 479 distributed across all phylostrata for both up-and down-regulated genes. Among them, the 480 oldest phylostratum (prokaryotic) showed the highest enrichment for the downregulation 481 response, being significantly enriched in >50% of the experiments. Furthermore, the 482 Gymnospermae phylostratum, representing gene families that appeared in the ancestor of seed 483 plants, tends to be upregulated in heat, salt, low iron, and phosphate. Conversely, species-484 specific genes showed the least enrichment, suggesting that newly formed genes are less 485 involved in the stress response. We conclude that stress responses tend to be regulated 486 independently from the age of the gene families, except the oldest (prokaryotic, more 487 responsive than expected) and the youngest (species-specific, less responsive than expected). 488 489

Discussion: 490
Cyanophora paradoxa is the only glaucophyte representative with a sequenced genome, and 491 due to the importance of this species for evolutionary studies we generated a comprehensive 492 expression atlas comprising several abiotic stresses. We focused the atlas on a wide range of 493 abiotic stresses for two reasons. The first is that co-expression networks require gene 494 expression data capturing as many different observations as possible (Usadel et al., 2009). The 495 second is that all organisms have developed transcriptional programs to cope with various 496 abiotic stresses (Zhu, 2016), but currently no study is available on gene expression comparison 497 at plant kingdom-wide scale. 498 The Cyanophora gene expression atlas revealed that there was no correlation between 499 the severity of growth inhibition caused by the stress and the number of transcriptional 500 changes (given as the number of differentially expressed genes) (r=0.278, p-value=0.437, Fig.  501 1e), even though both culture growth and the transcriptome were affected by the applied 502 stresses. For instance, stresses that cause a mild growth retardation (e.g. nitrogen and sulphate 503 deprivation) showed a high number of DEGs, while stresses causing a severe growth 504 inhibition (salt and mannitol) showed the lowest amount of DEGs. The overall low correlation 505 between growth inhibition and transcriptome changes suggests highly diverse strategies and 506 regulatory programs that control the stress acclimation. 507 We expanded the CoNekT-Plants database with the newly generated expression dataset 508 to give easy access to study gene expression and co-expression networks through the use of 509 various comparative tools. These tools allow us to overcome the paucity of knowledge related 510 to gene function in C. paradoxa. We exemplify these tools by studying phycobilisome 511 formation (Fig. 2c) and high light response (Fig. 3c), demonstrating how the co-expression 512 network can be used to predict the function of unknown genes. Furthermore, the updated 513

CoNekT-Plants can be used for cross-species analyses, as we have demonstrated by 514
identifying a conserved module involved in cell division in C. paradoxa and C. reinhardtii 515 (Fig. 4c). These tools thus allow to uncover the similarities and differences in gene expression 516 and functional modules across evolutionary distances, making CoNekT-Plants a unique 517 resource to study the evolution of functional modules. 518 Moreover, we analyzed whether the responses to abiotic stresses are conserved across 519 conditions and organisms, and we observed a conserved downregulation response related to 520 genes involved in photosynthesis and primary metabolism (Fig. 5c, 6). This observation is also 521 supported by the analysis of the most conserved families operating upon stress (Fig. 6b), in 522 which we find several, mostly downregulated photosynthesis and primary metabolism-related 523 families. Similar trends have been found in grasses, where the various members of the 524 Pooideae family showed highly diverse transcriptomic responses to cold stress, but a common 525 downregulation of transcripts involved in photosynthesis and metabolism (Schubert et al., 526 2019). Analysis of photosynthesis parameters during salt stress revealed that this is reflected 527 by a decrease of the total performance index of photosynthesis in multiple species (Pavlović et 528 al., 2019), showing a correlation between the transcript responses and plant phenotype. Our 529 findings suggest that photosynthesis is downregulated to prevent photoinhibition and cellular 530 damage not only in cold and salt stress but for a wider variety of stresses. Furthermore, since 531 this response is seen in algae and land plants, we propose that downregulation of 532 photosynthesis and metabolism is a kingdom-wide, stress-responsive program. 533 In contrast to the conserved downregulation responses, we observed only a modest 534 conservation of upregulation responses to various stresses across different species (Fig. 6a). 535 This is also reported for salt stress in six Lotus accessions, where only 1% of genes showed a 536 conserved response (Sanchez et al., 2011), in two strawberry cultivars, which displayed a 537 modest conservation of DEGs to the same pathogen (Wang et al., 2017), in Poaceae, which 538 revealed a poor expression conservation of orthologs across same tissue types (Davidson et al., 539 2012), and in seven Arabidopsis accessions which showed a divergent response of genes to 540 exogenous salicylic acid (van Leeuwen et al., 2007). We propose two non-exclusive 541 hypotheses for this surprising lack of conservation in transcriptomic responses. First, assuming 542 that orthologous sets of genes are still needed for acclimation to a given stress, investigating 543 only one layer of the gene activity regulation (transcript levels) might be insufficient. Second, 544 a high rate of rewiring, divergence, and the presence of multiple signal transduction pathways 545 might confound the signal. To provide a deeper understanding of stress responses, we suggest 546 that snapshots of additional regulatory layers, preferably closer to the end-product of gene 547 regulation (i.e., the protein), should be included in future studies. 548 Finally, we demonstrated that the oldest gene families, in contrast to the youngest 549 species-specific families, are actively responding to abiotic stress (Fig. 7). This suggests that 550 the young gene families are not typically involved in stress acclimation, which is in line with 551 the fact that the investigated stresses, and the corresponding coping mechanisms, are ancient.  Table S1. Composition of the C media. 804 Table S2. LSTrAP mapping statistics. The first column represents the sample names, the 805 second, third and fourth columns represent the reads mapped to the genome, the number of 806 reads mapped to noncoding genes ('no feature') and not uniquely ('ambiguous'), respectively. 807 The fifth, sixth, and seventh columns indicate the relative percentages. The eighth column 808 indicates the percentage of reads mapped to coding genes. 809 The first column shows the gene name, the second column indicates the condition the gene 813 resulted differentially expressed, the third and fourth columns indicate the adjusted p-value 814 and the log2 fold change associated to the gene and the experiment. 815 Table S8. Cyanophora paradoxa TPM-normalized gene expression matrix. 816 Table S9. HCCA clusters of C. paradoxa 817 Table S10