Evolution of cold acclimation and its role in niche transition in the temperate grass subfamily Pooideae

1

Evolution of cold acclimation in Pooideae

1

Corresponding author: Siri Fjellheim, Faculty of Biosciences, Norwegian University of Life Sciences, PO 2

BOX 5003, N-1432 Aas, Norway. Tel: +47 67232801 3

E-mail: [email protected] 4

5

2

Evolution of cold acclimation and its role in niche transition in the

6

temperate grass subfamily Pooideae

7

Marian Schubert

, Lars Grønvold

, Simen R. Sandve

, Torgeir R. Hvidsten

and Siri Fjellheim

8

*

Contributed equally 9

1

Department of Plant Sciences, Norwegian University of Life Sciences, Ås NO-1432, Norway.

10

2

Faculty of Chemistry, Biotechnology and Food Science, Norwegian University of Life Sciences, NO-1432, 11

Ås, Norway.

12

3

Centre for Integrative Genetics (CIGENE), Department of Animal and Aquacultural Sciences, Norwegian 13

University of Life Sciences, NO-1432, Ås, Norway.

14

4

Umeå Plant Science Centre, Department of Plant Physiology, Umeå University, SE-90187, Umeå, 15

Sweden.

16

Short summary: We show that cold acclimation is common in Pooideae, but that this adaptation to 17

temperate climate largely has evolved independently in different tribes.

18

Author contributions: All authors designed the experiment. M.S. and S.F. performed the growth 19

experiments. M.S. sampled and prepared RNA for sequencing, helped designing the data analysis pipeline, 20

contributed to the positive selection analysis and performed the phylogenetic analyses. L.G. developed, 21

implemented and conducted the transcriptome data analyses. All authors interpreted the results and wrote 22

the manuscript.

23

Funding information: The research was funded by grants from the Nansen Foundation to S.F. and the 24

TVERRforsk program at the Norwegian University of Life Sciences (NMBU) to S.F., T.R.H. and S.R.S.

25

This work was part of the PhD projects of M.S. and L.G. funded by NMBU.

26

Corresponding author e-mail: [email protected] 27

Keywords: cold acclimation, temperate adaptations, Pooideae, grasses, transcriptome, gene family 28

evolution, phylogeny 29

30

3

Abstract

31

The grass subfamily Pooideae dominates the grass floras in temperate regions, and has evolved complex 32

physiological adaptations to cope with environments characterized by extreme conditions like frost, winter 33

and strong seasonality. One such adaptation is cold acclimation, where plants increase their frost tolerance 34

in response to gradually falling temperatures and shorter days in the autumn. However, understanding how 35

complex traits like cold acclimation evolve through adaptive changes remains a major challenge in 36

evolutionary biology.

37

Here, we investigate the evolution of cold acclimation in Pooideae, and test a phylogenetically diverse set 38

of Pooideae species for cold acclimation capacity. To establish if the underlying genetic basis is shared, we 39

assembled transcriptomes of five phylogenetically diverse species and compared differential expression 40

after short- and long-term cold treatment. Furthermore, we studied the association between evolution of 41

gene family size and selection pressure on coding sequence with gene cold responsiveness.

42

All species have cold acclimation capacity, but the underlying genetic responses reflect widespread species- 43

specific cold responses in genes with conserved sequence. We saw evidence of protein coding sequence 44

and regulatory evolution, as well as the origin of novel genes and functions playing a role in cold response 45

evolution in Pooideae.

46

Our results are consistent with phylogenetic dating and paleoclimatic data suggesting that selection pressure 47

resulting from dramatic global cooling must have acted on already diverged lineages. However, phylogeny- 48

wide conservation of certain genes indicate that the Pooideae ancestor may have possessed some molecular 49

machinery to mitigate cold stress.

50

Introduction

51

Frost is one of the most dramatic abiotic stresses a plant can experience and exposure to frost is limiting 52

the opportunity for many species to diversify in temperate and arctic regions. Only a few, ancestrally 53

tropical angiosperm lineages have colonized temperate biomes (Judd et al., 1994; Wiens and Donoghue 54

2004; Kerkhoff et al., 2014). These lineages represent several independent transitions from tropical to 55

temperate environments and represent a variety of evolutionary adaptations to both unpredictable, short 56

incidences of episodic frost and predictable periods of periodic frosts (i.e. winter). Today’s cold climate is 57

the result of gradual cooling over the last 50 million years (my), and the temperate biomes originated 58

relatively recently in Earth’s history when the global climate started to cool in the late Eocene (Potts and 59

Behrensmeyer, 1992; Morley, 2000; Zachos et al., 2001; Fine and Ree, 2006; Eldrett et al., 2009; Liu et al., 60

2009; Strömberg, 2011). Around the Eocene-Oligocene (E-O) transition 34 my ago, global temperatures 61

fell dramatically (Pound and Salzmann, 2017). Throughout the Oligocene the temperatures continued to 62

drop and together with increased seasonality (Zachos et al., 2001; Eldrett et al., 2009) this triggered the 63

expansion of temperate biomes to intensify. The successful colonizers of the emerging temperate climates 64

faced gradually more severe stresses and evolved adaptations to frost, increased temperature seasonality, 65

and short growing seasons (Zachos et al., 2001; Eldrett et al., 2009; Mudelsee et al., 2014). The fact that 66

only a restricted number of plant lineages have transitioned into the temperate region suggests that it is 67

challenging to evolve the coordinated set of physiological changes needed to withstand low temperatures 68

(Donoghue, 2008).

69

The capacity for surviving environments with sub-zero temperatures depends on acquisition of frost 70

tolerance. Most importantly, plants need to maintain the integrity of cell membranes to avoid osmotic stress 71

during prolonged freezing (Thomashow, 1999). To endure predictable, periodic frost, plants must have 72

evolved to integrate both phenology and complex physiological adjustments. Through a process called cold

73

4 acclimation, freezing tolerance is acquired through a range of physiological changes governed by diverse 74

molecular pathways. These changes result in an increase in the sugar content of cells, change in lipid 75

composition of membranes and synthesis of anti-freeze proteins (Janska et al., 2010; Preston and Sandve, 76

2013). In addition, low non-freezing temperatures may affect plant cells by decreasing metabolic turnover 77

rates, inhibiting the photosynthetic machinery and decreasing stability of biomolecules (e.g. lipid 78

membranes) (Sandve et al., 2011; Crosatti et al., 2013). Temperate and arctic plants use the gradually lower 79

temperature and day-length in the autumn as cues to initiate cold acclimation.

80

In the grass family (Poaceae), the subfamily Pooideae is the one occupying the coldest climate space 81

(Edwards and Smith, 2010). Except for a few hundred species in early diverging Pooideae tribes, most of 82

the ca. 4200 Pooideae species belong to the species-rich ‘core Pooideae’ clade (sensu Davis and Soreng, 83

1993; Soreng and Davis, 1998). The closest sister group to the core Pooideae is the tribe Brachypodieae, 84

which contains the model grass Brachypodium distachyon (Soreng et al., 2015). Because core Pooideae 85

contain economically important species like wheat (Triticum aestivum) and barley (Hordeum vulgare) as 86

well as several forage grasses like ryegrass (Lolium perenne) there is extensive knowledge about 87

physiological and molecular mechanisms underlying frost tolerance, however, only from a limited number 88

of core Pooideae species (e.g. Thomashow, 1999; Sandve et al., 2008; Galiba et al., 2009; Sandve et al., 89

2011, Preston and Sandve 2013; Fjellheim et al., 2014; McKeown et al., 2016; Woods et al., 2016;

90

McKeown et al., 2017). It has been shown in numerous studies that several species of the core Pooideae 91

have the capacity to cold acclimate. Facilitating these responses are, among others, five gene families 92

known to play important roles during cold-stress response and cold acclimation in core Pooideae. These 93

gene families code for C-repeat binding factors (CBF), dehydrins (DHN), chloroplast-targeted cold- 94

regulated proteins (ctCOR), ice recrystallization inhibition proteins (IRIP) and fructosyl transferases (FST).

95

(CBF: Badawi et al., 2007; Li et al., 2012. DHN: Olave-Concha et al., 2004; Rorat, 2006; Kosová et al., 96

2007; 2014. ctCOR: Gray et al., 1997; Crosatti et al., 1999; 2013; Tsvetanov et al., 2000. IRIP: Antikainen 97

and Griffith 1997; Hisano et al., 2004; Kumble et al., 2008; Sandve et al., 2008; John et al., 2009; Zhang et 98

al., 2010; Sandve et al., 2011. FST: Hisano et al., 2004; Tamura et al., 2014).

99

Although Pooideae species dominate temperate and arctic grass floras (Hartley, 1973; Visser et al., 2014), 100

the ancestors of this group were most likely adapted to tropical or subtropical climates (Bouchenak-Khelladi 101

et al., 2010; Strömberg 2011). Pooideae originated in the late Cretaceous or early Paleogene period 102

(Bouchenak-Khelladi et al., 2010; Christin et al., 2014; Spriggs et al., 2014; Schubert et al., 2018) during a 103

time where the global climate in general was warm (Zachos et al., 2001; Mudelsee et al., 2014) and 104

seasonality in temperature relatively low (Archibald et al., 2010). The colonization of temperate biomes by 105

Pooideae was likely facilitated by evolution of successful adaptations to cold conditions. Although studies 106

have inferred adaptation to cooler environments at the base of the Pooideae phylogeny (Edwards and Smith, 107

2010; Schubert et al., 2018), it is still not known whether the Pooideae’s most recent common ancestor 108

already was adapted to cold and/or freezing conditions, or if adaptations to cold evolved during 109

diversification of Pooideae along with falling temperatures throughout the Cenozoic and expansion of 110

temperate biomes.

111

In this paper, we aim to investigate the evolution of frost tolerance in response to cold acclimation in the 112

grass subfamily Pooideae. We first conduct classical cold acclimation and freezing test experiments to 113

determine if cold acclimation is a shared trait throughout the Pooideae phylogeny. Next, we take advantage 114

of comparative transcriptomics analyses to investigate if a shared molecular basis of cold acclimation in 115

Pooideae is present. Finally, we interpret our results in a phylogenetic and paleoclimatic context, and 116

suggest a model of mostly independent evolution of cold acclimation responses.

117

5

Results

118

Cold acclimation exists in early diverging lineages

119

We conducted a classical cold acclimation and freezing test experiment in nine species that represent major, 120

species rich clades in the Pooideae subfamily or belong to very early diverging lineages (Fig. 1A), and that 121

are distributed in areas where they are regularly exposed to cold, sub-zero temperatures (Fig. 1B). Freezing 122

tests revealed that cold acclimation, i.e. increased frost tolerance through exposure to cold non-freezing 123

periods, exists both in core Pooideae and species of early diverging lineages (Fig. 2). All acclimated plants 124

exhibited higher regrowth capacity at -4 and -8°C compared to non-acclimated plants, although the increase 125

in regrowth capacity was not significant at -4°C for S. lagascae and N. stricta. Non-acclimated plants of 126

early diverging Pooideae species performed better at -4°C than non-acclimated Brachypodium species and 127

H. vulgare, and comparable with the regrowth capacity of non-acclimated plants of the perennial core 128

Pooideae species L. perenne and Elymus repens.

129

De novo transcriptome assembly identified 8633 high confidence ortholog groups

130

To investigate the molecular basis of cold response evolution, we sampled leaf material from five species 131

spanning early to later diverging Pooideae lineages before and after subjecting them to a drop in temperature 132

and shorter days (Fig. S1). The transcriptome of each species was assembled de novo from RNA-sequencing 133

(RNA-Seq) reads (Table S1). Ortholog groups (OGs) were inferred by using protein sequences from the 134

five de novo assemblies and the reference genomes of L. perenne, H. vulgare, B. distachyon, Oryza sativa, 135

Sorghum bicolor and Zea mays. A set of 8633 high confidence ortholog groups (HCOGs) was identified 136

after filtering based on gene tree topology and species representation (Table S2).

137

Shared cold response genes included known abiotic stress genes

138

To investigate shared expression patterns across species, we constructed a single cross-species gene 139

expression table with HCOGs as rows and samples as columns - by summing the expression values of 140

monophyletic paralogs and setting the expression of missing orthologs to zero (Table S3). After removing 141

differences in mean expression levels between species for each gene, clustering reconstructed a tree with 142

replicates, and then time points, grouping together (Fig. S2A). The fact that time points mostly clustered 143

before species indicated a common transcriptional response to cold across species. An exception was time 144

points W4 and W9, which tended to cluster together and by species, indicating that responses after 4 and 9 145

weeks were very similar. We also observed a clear effect of the diurnal rhythm, with time points sampled 146

in the morning (W0, W4 and W9) forming one cluster and time points sampled in the afternoon (D0 and 147

D1) forming another.

148

Next, we focused on the short- and long-term cold response by analyzing changes in gene expression from 149

before cold treatment to eight hours (short-term) and 4-9 weeks (long-term) of cold treatment. For all 150

species pairs, there was a low, but statistically significant, correlation between expression changes after 151

cold treatment of orthologs in HCOGs (Fig. S2B). This indicates that orthologs tending to change 152

expression in the same direction after cold treatment, an observation that explains the clustering of samples 153

according to treatment before species (Fig. S2A). Finally, we analyzed the statistical significance of 154

expression change after cold treatment by classifying each individual gene as either differentially expressed 155

or not (FDR adjusted p-value < 0.05 and fold change > 2 or < 0.5, Table S4). Although a considerable 156

number of genes responded to cold in each species (1000-3000 genes), the overlap between species were 157

low. Of the 5577 HCOGs with at least one differentially expressed gene (DEG), 50% contained only one 158

DEG and thus represented orthologs differentially expressed in only one species. 31% contained DEGs 159

from two species, 14% from three species and only 4% contained DEGs from four species. Importantly, 160

this pattern was independent of expression level and did not change for the subsets of HCOGs with medium

161

6 or high expression (Table S4). Sixteen genes shared the same cold response (either short- or long-term) in 162

the same direction (either up or down) in all five species (Table 1). These genes thus represent strong 163

candidates for responses to cold that might have been conserved throughout the evolution of Pooideae. Nine 164

of these genes belonged to families known to be involved in cold stress or other abiotic stress responses in 165

other plant species. Of the sixteen conserved genes, 12 genes displayed short-term up regulation indicating 166

that stress response, as opposed to long-term acclimation response, is potentially more conserved.

167

We compared our DEGs to a compilation of 55 H. vulgare genes shown to be responsive to low temperature 168

in several previous microarray experiments (Greenup et al., 2011). We could map 33 of these genes to 169

unique OGs, of which 11 were HCOGs. We observed significant similarity in cold response between the 170

33 previously identified cold response genes and the short-term DEGs in our data (P < 0.05, Fig. S3). This 171

similarity was statistically significant for all five species, although noticeably larger in H. vulgare than in 172

the other species. This comparison thus shows that our transcriptome data was consistent with previous 173

findings in H. vulgare, and that cold response genes identified in H. vulgare exhibit some cold response in 174

other Pooideae.

175

Biological processes involved in cold response

176

To identify biological processes that evolved regulation at different stages of Pooideae speciation, we 177

targeted gene sets that were exclusively differentially expressed in all species within a clade in the 178

phylogenetic tree (i.e. branch specific DEGs), and tested these for enrichment of Gene Ontology (GO) 179

biological process annotations (Fig. 3A). For the genes that were differentially expressed in all our species 180

(Pooideae base [PB]), we found that up-regulated genes were enriched for GO annotations related to signal 181

transduction (two pseudo response regulators and diacylglycerol kinase 2 (DGK2)) and abiotic stimulus 182

(Gigantea, LEA-14, DnaJ and DGK2), while down-regulated genes were enriched for photosynthesis and 183

metabolism. For the genes that were exclusively differentially expressed in all species except N. stricta 184

(early split [ES]), down-regulated genes were again enriched for metabolism and photosynthesis.

185

Positive selection in cold responsive genes

186

To resolve the uncertain placement of Meliceae and Stipeae lineages within the Pooideae species tree we 187

reconstructed gene trees using the 3914 HCOGs with exactly one sequence from each of the five Pooideae 188

species and rice. In the most common gene tree topology, S. lagascae and M. nutans formed a monophyletic 189

clade (755 trees), but topologies where either M. nutans or S. lagascae diverged first (434 and 271 trees, 190

respectively) were also common (Fig. S4). Averaging the posterior probabilities of the branches leading to 191

the MRCA of S. lagascae and M. nutans did not produce significant support for any of the four most 192

common topologies. Thus, we displayed the relationships of Meliceae and Stipeae tribes as unresolved 193

trichotomy in Fig. 3 and 4, which is congruent with 1460 of the gene trees.

194

We also tested for positive selection (i.e. selection for beneficial substitutions) in coding sequences at each 195

of the internal branches of the species tree. The tests were only performed on the branches where the 196

respective HCOG gene tree topology was compatible with the species tree topology. 16-18% of the HCOGs 197

showed significant signs of positive selection (P < 0.05) depending on the branch (Fig. 3B). Next, we tested 198

for overrepresentation of positive selection among the branch specific DEGs. Although not statistically 199

significant, there was a tendency that gain of cold response was associated with positive selection on coding 200

sequence at the early (ES) and late split (LS) branches (P = 0.077 and P = 0.072, respectively) (Fig. 3B).

201

This observation indicates that both selection on regulatory as well as protein coding sequence may have 202

played a role in evolving cold responses in Pooideae.

203

204

7

Cold acclimation gene families have diverse evolutionary histories

205

Of the 33 previously described cold responsive H. vulgare genes that were also differentially expressed in 206

our experiment (Fig. S3), as many as 22 were not included in the HCOGs, mainly because they belonged 207

to gene families with complex gene duplication histories in lineages leading up to two or more of our studied 208

species. Next, we tested if this was a general trend for all cold responsive genes. This analysis showed that 209

cold responsive genes in general are members of larger gene families than genes that do not respond to 210

cold, with the exception of short-term down-regulated genes (Fig. S5). Conversely, short-term down- 211

regulated genes were enriched in single-copy gene families while up-regulated genes were less common in 212

these families (Fig. S5). We also tested if response to cold was associated with the age of the gene family 213

through a phylostratigraphic analysis. The general trend was that cold-responsive genes tended to be 214

underrepresented among the most ancient families (Fig. S6).

215

We conducted detailed analyses of five of the best studied cold acclimation gene families coding for C- 216

repeat binding factors (CBF), dehydrins (DHN), chloroplast-targeted cold-regulated proteins (ctCOR), ice 217

recrystallization inhibition proteins (IRIP) and fructosyl transferases (FST). We reconstructed gene trees 218

for these families (Fig. S7-S14) using reference sequences and de novo assembled transcripts (Table S5), 219

and overlaid information on short- and long-term cold response DEGs. Dehydrin genes are well studied in 220

H. vulgare with 13 known dehydrin homologs (termed HvDHN1 to HvDHN13). Structurally they can be 221

grouped into four distinct types based on the presence of amino acid segments (Y, K and S): SK

-type 222

(HvDHN8), KS-type (HvDHN13), K

-type (HvDHN5) and Y

SK

-type dehydrins (the 10 remaining 223

HvDHN homologs) (Kosová et al. 2007). Because the genes for these four groups represent 224

phylogenetically distinct clades (Karami et al. 2013), we reconstructed individual gene trees for each group, 225

however not for HvDHN5, because homologs were only found in the H. vulgare transcriptome. Two groups 226

of CBF genes – CBFIIIc/d and CBFIV – are especially important for cold acclimation in Pooideae and are 227

restricted to this subfamily (Badawi et al. 2007, Li et al. 2012). As the CBF gene family is large and highly 228

complex, we restricted our analyses of the CBF gene family to these two groups. Our analyses could not 229

completely resolve the topologies for most of the gene trees owing to their complex evolutionary history.

230

We identified, however, that DHN8, DHN13 of the DNH gene family and the ctCOR gene family as likely 231

candidates for conserved cold acclimation genes within the Pooideae (Fig. 4). Also, the expansion of the 232

CBFIIId group likely started in the earliest Pooideae lineage. In contrast, the evolution of the IRI and FST 233

motifs, the expansion of the CBFIV group and the evolution of DHN5 (no tree reconstructed) was restricted 234

to later diverging lineages. Furthermore, these detailed analyses revealed a decreasing similarity in 235

transcriptional cold response with increasing phylogenetic distance (Fig. 4). Naturally, this analysis is 236

biased towards H. vulgare as most of what we know about these gene families is derived from crop species.

237

Discussion

238

Cold acclimation evolved independently in separate lineages

239

It has been suggested that the ecological success of the Pooideae subfamily in the northern temperate 240

regions critically relies on adaptations to cold temperatures. Here we show that cold acclimation increased 241

frost tolerance relative to non-acclimated plants in all the major Pooideae lineages (Fig. 2). This 242

demonstrates that the cold acclimation trait is pervasive across Pooideae, and an important part of 243

adaptation to a life in temperate and arctic biomes.

244

To dissect the molecular basis of the cold acclimation responses and determine if the cold acclimation 245

response has a shared molecular basis, we conducted a global comparative, transcriptomic analysis for five 246

phylogenetically diverse species. For each species, 10-30% of the highly conserved genes (i.e. genes in the 247

8,633 HCOGs) responded to cold, which is in line with what we know about cold acclimation programs

248

8 from other species (e.g. Arabidopsis (Park et al., 2015)). Transcripts involved in photosynthesis and 249

response to abiotic stimuli were significantly enriched among the genes with cold response in all species 250

(Fig. 3A). It is known that down-regulation of the photosynthetic machinery during cold temperatures is 251

one mechanism to prevent photoinhibition (Crosatti et al., 2013) and cellular damage. This protective 252

mechanism might therefore have existed early in the evolution of Pooideae. Enrichment of genes involved 253

in photosynthesis among differentially expressed genes is not specific to cold stress, but can be found as 254

response to other stresses as well, among them salt stress (Yamamoto et al., 2015) as well as heat and 255

drought stress (Liu et al., 2015). Interestingly, in our study nearly half of all cold responsive genes were 256

species-specific in differential expression tests (Table S4). This pattern was also supported by low (although 257

significant) correlations in gene expression fold-changes (Fig. S2B). These results indicate that evolution 258

of cold acclimation evolved at least partly independently in separate lineages and are in line with recent 259

work on transcriptional responses to vernalization in Pooideae (McKeown et al., 2016; Zhong et al., 2018).

260

Many genes previously associated with cold acclimation belong to complex gene families (Thomashow, 261

1999; Sandve and Fjellheim, 2010) and are therefore not included in the global analysis of highly conserved 262

orthologs. We therefore conducted an in-depth analyses of gene families known to be important for cold 263

acclimation (Fig. 4, Fig. S7-14) in well studied core Pooideae species (Choi et al., 2002; Vágújfalvi et al., 264

2003; Badawi et al., 2007; Knox et al., 2008; Zhang et al., 2010; Livingston et al., 2009; Knox et al., 2010;

265

Soltesz et al., 2013; Jeknić et al., 2014; Todorovska et al., 2014; Marozsán-Tóth et al., 2015). The lack of 266

conserved expression patterns of these genes across the tested Pooideae species (Fig. 4) further supports 267

that cold acclimation is regulated differently in the five Pooideae species.

268

Three gene families or gene family members known to be involved in cold acclimation in core Pooideae 269

crop species evolved de novo or evolved a new function in some lineages. Firstly, the FST gene family 270

originated at the base of the core Pooideae (Fig. S14). Secondly, the functional domain of the IRIP-gene 271

family evolved at the base of core Pooideae/Brachypodium (Fig. S13). Lastly, in the DHN gene family the 272

gene DHN5 evolved in Triticeae (Fig. 4). Together, these three examples provide conclusive evidence that 273

the cold acclimation responses in early diverging lineages at least partly rely on other genes and mechanisms 274

than the core Pooideae and Brachypodium. Some very likely candidates for cold responses specific to early 275

diverging Pooideae lineages are homologs of the Y

SK

-type DHN genes in the DHN family (Fig. 4, Fig.

276

S11). Y

SK

-type DHN genes are neither known to be involved in cold acclimation of Brachypodium and 277

core Pooideae, nor have we found any support for that in our data.

278

Even though our analyses show that the species largely follow independent evolutionary trajectories, the 279

responses to a large extent build on a common genetic basis. The CBF, DHN and COR gene families are 280

known to be involved in cold stress adaptation across the plant kingdom (Thomashow 2001; Rorat 2006) 281

and FST genes has evolved several times from vacuolar invertases independently across the angiosperms 282

(Vijn and Smeekens 1999). This indicate that some genes code for proteins with biochemical functions 283

suited to be recruited for cold stress mitigation. One gene family whose function in cold acclimation might 284

be conserved throughout the Pooideae is the ctCOR gene family (Fig. 4, Fig. S12), rendering this gene 285

family an interesting candidate that could constitute a shared molecular basis of the cold acclimation trait 286

in Pooideae. Interestingly, a ctCOR homolog is upregulated as response to short-term cold in rice 287

(Maruyama et al., 2014), indicating that the co-option of this gene to the highly specialized cold acclimation 288

pathway in Pooideae build on a stress response to cold that is conserved well beyond Pooideae. A very 289

similar scenario is likely for the two genes DHN8 and DHN13. Homologs of both genes are induced during 290

chilling stress in rice (Aguan et al., 1991; Lee et al., 2005) and there is strong support for both genes having 291

a function in cold acclimation of core Pooideae (Kosová et al., 2007). Our results suggest that DHN8 and 292

DHN13 are single copy genes involved in cold acclimation in all but one investigated species.

293

9 Although the molecular mechanisms behind adaptive evolution is still an active field of research, it is now 294

indisputably established that novel gene regulation plays a crucial role (Romero et al., 2012). The evolution 295

of gene regulation proceeds by altering non-coding regulatory sequences in the genome, such as (cis-) 296

regulatory element (Wittkopp et al., 2012), or by altering the coding sequences of regulatory proteins. The 297

high number of genes with species-specific cold responses observed in this study is consistent with the 298

recruitment of genes with existing cold tolerance functions by means of regulatory evolution. However, 299

previous studies have also pointed to the evolution of coding sequences (Vigeland et al., 2013) as underlying 300

the acquisition of cold tolerance in Pooideae. To investigate possible adaptive evolution in coding regions, 301

we tested for the enrichment of positive selection among branch specific cold responsive genes (Fig. 3B).

302

Although not statistically significant (P = 0.07), increased signals for positive selection on coding 303

sequences was observed for cold responsive genes in a period of gradual cooling preceding the E-O event.

304

In conclusion, we find evidence to support that increased selection for protein coding sequence, changes in 305

gene regulation as well as novel genes and functions play a role in evolution of cold responses in Pooideae.

306

In general, it is important to note that lack of shared signals between lineages can also arise for other reasons 307

than biological meaningful differences and must therefore be interpreted with caution. In our study, limited 308

statistical power and the use of de novo transcriptomes could impact our ability to detect overlap in cold 309

responses between species. To minimize the effect of these putative biases in this study we only analysed 310

ortholog groups containing genes that could be assembled across all species in the phylogeny, ensuring that 311

we measured changes in gene expression levels and not changes in gene content. Furthermore, since three 312

of the species lacked reference genomes, we employed a de novo assembly pipeline to reconstruct the 313

transcriptomes. We showed that this pipeline could recover a set of H. vulgare genes previously identified 314

as cold responsive (Greenup et al., 2011) and that most of these genes were also cold responsive in H.

315

vulgare in our experiment (Greenup et al., 2011). Moreover, we conducted analyses to assess if sequencing 316

depth was a significant factor in our ability to identify shared patterns of differential expression (Table S4).

317

These results clearly show that sequencing coverage cannot explain the high ratio of lineage-specific cold 318

acclimation responses, indicating that increasing the statistical power will not change the conclusions of 319

this study. Hence, the global trend of low levels of shared cold acclimation responses is unlikely to be 320

methodological artefact, and consistent with the evolution of cold tolerance mechanisms being shaped 321

mostly independently in the different Pooideae lineages.

322

A shared adaptive potential for evolution of cold adaptation

323

Our differential expression analysis revealed 16 genes sharing cold associated gene expression shifts across 324

all species (Table 1). An overwhelming majority (12/16) of the conserved genes were short-term cold 325

responsive (Table 1) and this observation suggests that existing stress genes have been the first to be co- 326

opted into the cold response program. Ten of these genes have previously been shown to be regulated by 327

other abiotic and biotic stressors (Table 1). Nine of these conserved cold response genes are involved in 328

response to abiotic stresses in other plants, such as osmotic stress and drought. The SNAC1-like/NAC 329

transcription factor 67 is one example, with homologs induced by osmotic stress in the three monocots 330

Musa × paradisiaca (Negi et al., 2015; Yang et al., 2015), Oryza sativa (Nakashima et al., 2012) and 331

Triticum aestivum (Mao et al., 2014). Co-option of such genes into a cold-responsive pathway might have 332

been the key to acquire cold tolerance. Some responses may be even more ancient than the MRCA of the 333

Pooideae, as some of the genes with cold response in all five species are also expressed in response to cold 334

in other species (Table 1), amongst them rice, which is a related tropical species with some tolerance to 335

chilling stress (Wang et al., 2016). Shaar-Moshe et al., (2017) recently identified B. distachyon genes 336

induced by drought, high salinity or heat as well as by their combinatorial effects. All genes from Table 1 337

were differentially expressed in at least one of their stress conditions. Noticeably however, none of those 338

genes was responsive to drought stress alone. We also found that cold responsive genes are

339

10 underrepresented among the most ancient gene families (Fig. S6). Taken together, it seems as the 340

transcriptional cold responses that are shared across the Pooideae phylogeny represent deeply conserved 341

ancestral responses to mitigate general cellular stress, rather than representing specific adaptations to 342

prepare for periodic frost.

343

Pooideae lineages evolved specific cold adaptation by expanding gene families

344

Gene family expansion has previously been suggested to play a role in cold adaptation in Pooideae (Sandve 345

and Fjellheim 2010; Li et al., 2012). As previously discussed, the conservative filtering of ortholog groups 346

employed in this study removed complex gene families containing duplication events shared by two or 347

more species. A majority of the previously described H. vulgare cold responsive genes (Greenup et al., 348

2011) belonged to gene families with complex gene duplications histories. Moreover, an analysis of all cold 349

responsive genes confirmed a positive association between gene family size and cold response in our data, 350

with the excepttion of short-term down-regulated genes which instead were enriched among single copy 351

gene families (Fig. S5). These results thus confirm previous findings that stress related genes are enriched 352

in fast evolving gene families (Panchy et al., 2016). The gene trees of the gene families targeted for 353

phylogenetic analyses (Fig. 4, Fig. S7-14) displayed several duplication events and gene family expansions.

354

Based on our results we propose that gene family expansion was an important mode of cold adaptation in 355

the Pooideae subfamily. Our results corroborate findings from Sandve and Fjellheim (2010), who identified 356

an increase of gene copy number in the CBFIV group and the FST and IRIP gene families as an evolutionary 357

force of cold climate adaptation of core Pooideae and Brachypodium species. Although we lack sufficient 358

genomic data for early diverging Pooideae species, we found evidence for Pooideae specific expansions in 359

the CBFIIId group the ctCOR gene family and in Y

SK

-type DHN genes. Expansion of gene families may 360

have led to functional specialization or novel functions of the various gene copies. Other studies have shown 361

that stress related gene families tend to expand via tandem duplications (Hanada et al., 2008), which may 362

lead to lineage-specific expansion of the gene family (Lespinet et al., 2002). Although de novo assembly of 363

transcriptomes from short-read RNA-Seq data is a powerful tool that has vastly expanded the number of 364

target species for conducting transcriptomic analysis, the approach has limited power to distinguish highly 365

similar transcripts such as paralogs. Homologous genes might exist in early diverging Pooideae species, 366

but may not be expressed. Furthermore, coding sequences may contain in-sufficient informative 367

substitutions to reconstruct the true species topology. Further insight into the role of gene duplication events 368

in evolution of the Pooideae cold adaptation would therefore benefit immensely from additional reference 369

genomes.

370

Evolution of cold adaptation in a paleoclimatic context

371

At face value, the most striking result from our study is the presence of cold acclimation capacity in all our 372

study species, and at the same time a generally low similarity between the species in transcriptomic 373

responses to cold acclimation, as well as the presence of well described cold acclimation genes, such as 374

FSTs, DHN5 and functional IRIP, in only some of the species. However, in a phylogenetic and 375

paleoclimatic context, these contrasting patterns can easily be reconciled. A recent study by Schubert et al.

376

(2018) places the divergence of all study species well before the E-O transition, when the global climate 377

was warm, seasonality low and temperate biomes had yet to expand. Thus, the most likely scenario is that 378

our study species have experienced >45 million years of independent evolution, and that adaptations to 379

periodic frost have independent evolutionary origins.

380

Responding to episodic and periodic frost are two fundamentally different physiological processes (Körner, 381

2016). Episodic frost induces a stress response and such existing responses could be exploited to tackle 382

frost when first experienced by the plants in the cooling climate during Eocene. On the other hand, 383

responses to periodic frost require controlled changes in cellular components initiated by plants upon

384

11 sensing shortening photoperiod and lower temperatures in the autumn. This complex response must have 385

evolved over a long period of evolutionary time, during which the plants were exposed to increasingly more 386

severe winters, exerting selection pressure for the cold acclimation response. Ancestral state reconstructions 387

of climatic niches of Pooideae species suggest that severe winter was not encountered by the Pooideae 388

ancestor, but was first experienced by Pooideae species more recently in the evolutionary history (Schubert 389

et al., 2018). Based on phylogenetic and paleoclimatic data we would therefore expect that selection 390

pressure for adaptations to winter seasons and frost spells across the Pooideae phylogeny would result in 391

independent evolution of cold acclimation capacity and many lineage and species specific gene expression 392

responses, as shown by our analyses.

393

Nevertheless, ancestral state reconstruction of the climatic niche of the ancestor shows that it likely 394

experienced episodic incidents of frost and that Pooideae possibly originated in a temperate micro-niche in 395

emerging alpine orogenies (Schubert et al., 2018). The conserved set of stress response genes identified in 396

our comparative transcriptome analysis may have represented a first tool-kit to cope with frost stress. This 397

may have represented a fitness advantage for the Pooideae ancestor in the newly emerging environment 398

with incidents of mild frost, allowing time to evolve the more complex physiological adaptations required 399

to endure the temperate climate with strong seasonality and cold winters that emerged following the E-O 400

transition (Eldrett et al,. 2009). Thus, preadaptations to temperate climate may have given Pooideae species 401

the advantage to expand and diversify with the expanding temperate biomes after the E-O transition, while 402

gradually building more complex adaptations to cold. This is in line with Schubert et al. (2018) that show 403

that diversification rates in Pooideae is inversely correlated with temperature.

404

Materials and Methods

405

Plant material

406

Seeds from nine Pooideae species were acquired from germplasm collections or collected in nature. The 407

selected species represent major, species rich clades in the Pooideae subfamily, or belong to very early 408

diverging lineages (Soreng et al., 2015). All species are distributed in areas where they are regularly 409

exposed to cold, sub-zero temperatures. Furthermore, to ensure that the accessions of the wild species that 410

were used in the experiment are adapted to cold conditions, accessions were selected from localities with 411

cold winters (mean temperature of the coldest quarter is below zero degrees, data extracted from the 412

WorldClim v.1.4 dataset (Hijmans et al., 2005) (Fig. S15A).

413

Nardus stricta L. (2n=2x=26) is a perennial species, distributed in Europe, western parts of Asia, and North 414

Africa and introduced to New Zealand and North America (Tutin, 1980; Hultén and Fries, 1986; Clayton 415

et al., 2006). Seeds were collected in July 2012 in Romania, [46.69098 N, 22.58302 E]. The locality of the 416

collected accession experiences sub-zero temperatures in November - March and the mean temperature of 417

the coldest quarter is -3°C (Fig. S15B).

418

Stipa lagascae Roem and Schult. (2n=2x=22) is a perennial species that is distributed in temperate regions 419

around the Mediterranean Sea and parts of temperate West Asia (Tutin, 1980; Hultén and Fries, 1986;

420

Clayton et al., 2006). Seeds from accession PI 250751 were acquired from U.S. National Plant Germplasm 421

System (U.S.-NPGS) via Germplasm Resources Information Network (GRIN). The accession was collected 422

49 km northwest of Tabriz, Iran in 1958. The locality of the collected accession experiences sub-zero 423

temperatures in December - March and the mean temperature of the coldest quarter is –0.2°C (Fig. S15B).

424

Melica nutans L. (2n=2x=18) is a perennial species distributed in temperate parts of Eurasia (Tutin, 1980;

425

Hultén and Fries, 1986; Clayton et al., 2006). Seeds were collected in Germany, [50.70708 N, 11.23838 E],

426

12 in June 2012. The locality of the collected accession experiences sub-zero temperatures in December - 427

March and the mean temperature of the coldest quarter is –0.9°C (Fig. S15B).

428

Brachypodium distachyon (L.) P. Beauv. (2n=2x=10) is an annual species natively distributed in Europe, 429

East Africa and temperate parts of West Asia (Tutin, 1980; Hultén and Fries, 1986; Clayton et al., 2006).

430

Seeds for accession ‘Bd1-1’ (W6 46201) were acquired from U.S.-NPGS via GRIN. This accession is a 431

singel-seed descent inbred line of the accession PI170218 (U.S.-NPGS, GRIN) collected in Soma, Manisa, 432

Turkey in 1948. The line is classified as a winter variety, i.e. Bd1-1 requires vernalization to flower and it 433

is freeze tolerant (Colton-Gagnon et al., 2014).

434

B. pinnatum (L.) P. Beauv. (2n=28) is a perennial species distributed in temperate parts of Eurasia (Tutin, 435

1980; Hultén and Fries, 1986; Clayton et al., 2006). Seeds were collected in October 2015, Norway 436

[59.71861 N, 10.59333 E]. The locality of the collected material experience sub-zero temperatures in 437

November - March, and the mean temperature of the coldest quarter is –3.0°C (Fig. S15B).

438

B. sylvaticum (Huds.) P. Beauv. (2n=2x=18) is a perennial species distributed in temperate parts of Eurasia 439

(Tutin, 1980; Hultén and Fries, 1986; Clayton et al., 2006). Seeds were collected in October 2015, Norway 440

[59.68697 N, 10.61012 E]. The locality of the collected material experience sub-zero temperatures in 441

November - March, and the mean temperature of the coldest quarter is –3.2°C (Fig. S15B).

442

Hordeum vulgare L. (2n=2x=14) seeds for cultivar ‘Igri’ were provided by Prof. Åsmund Bjørnstad, 443

Department of Plant Sciences, Norwegian University of Life Sciences, Norway. Igri is a two-row winter 444

barley developed in Germany in 1976 and is much used across Western Europe (Fischbeck, 2003). Igri 445

requires vernalization to flower and is tolerant to frost (Kosová et al., 2010). Sonja is a two-row winter 446

barley developed in Germany in 1974 (Fischbeck, 2003).

447

Lolium perenne L. seed for cultivar ‘Fagerlin’ were provided by Dr. Kovi, Department of Plant Sciences, 448

Norwegian University of Life Sciences, Norway. Fagerlin is a winter-hardy cultivar developed in Norway 449

in 2008.

450

Elymus repens (L.) Gould. (2n=4x 6x=28 42) is a perennial species distributed in Europe and Northern 451

parts of Asia and is introduced to North America (Tutin, 1980; Hultén and Fries, 1986; Clayton et al., 2006).

452

Seeds were collected in October 2015, Norway [59.66111 N, 10.89194 E]. The locality of the collected 453

material experience sub-zero temperatures in November - March, and the mean temperature of the coldest 454

quarter is –3.1°C (Fig. S15B).

455

Freezing tests

456

Seeds from all nine species (see Plant material, for Hordeum vulgare the ‘Sonja’ cultivar was used) were 457

germinated and plants grown in the green house at 20°C under natural day light. Each individual was 458

divided into four clones, one for each treatment and control. The plants were acclimated at 4°C and short 459

(8h) days for three weeks. Control conditions were short days and 20°C. The light intensity was 50 µmol 460

m

s

. At the end of the cold acclimation period, plants were subjected to freezing at three different 461

temperatures (-4, -8 and -12°C) following Alm et al., (2011). For each temperature we used 15 acclimated 462

and 15 non-acclimated individuals per species. After freezing, plants were cut down to approximately 3 cm 463

and grown at 20°C under long days in a greenhouse with natural light conditions. Two and three weeks 464

after the plants were moved into 20°C and long days they were assessed for regeneration ability and scored 465

from 0 (dead) to 9 (growth without damage). Differences between acclimated and non-acclimated 466

individuals within each species were tested with a one-tailed t test in R (R Core Team 2016) using the ‘stats’

467

package.

468

13

Cold treatment: Sampling, RNA extraction and sequencing

469

Seeds from N. stricta, M. nutans, S. lagascae, B. distachyon and Hordeum vulgare (cultivar ‘Igri’) were 470

germinated and initially grown in a greenhouse at a neutral day length (12 hours of light), 17°C and a 471

minimum artificial light intensity of 150 µmol/m

s. To ensure that individual plants were at comparable 472

developmental stages at sampling, the onset of treatment for different species was based on developmental 473

stages rather than absolute time. Most importantly, none of the plants had transitioned from vegetative to 474

generative phase, as this could have affected the cold response. Plants were grown until three to four leaves 475

had emerged for M. nutans, S. lagascae, B. distachyon and H. vulgare, or six to seven leaves for N. stricta 476

(which is a cushion forming grass that produces many small leaves compared to its overall plant size).

477

Depending on the species, this process took one (H. vulgare), three (B. distachyon and S. lagascae), six (M.

478

nutans) or eight (N. stricta) weeks from the time of sowing. Subsequently, plants from all species were 479

randomized and distributed to two cold chambers with short day (8 hours of light), constantly 6°C and a 480

light intensity of 50 µmol/m

s. Plants were kept in cold treatment for the duration of the experiment. Leaf 481

material for RNA isolation was collected i) in the afternoon (at zeitgeber (ZT) 8) on the day before cold 482

treatment (D0) and in the afternoon (ZT 8) on the first day of cold treatment (8 hours after initiation of cold 483

treatment, D1) and ii) in the morning (ZT 0) before cold treatment (W0), 4 weeks of cold treatment (W4) 484

and 9 weeks of cold treatment (W9) (Fig. S1). The sampling time points were chosen to be able to separate 485

chilling stress responses (first day of treatment) and long-term responses that represent acclimation to 486

freezing temperatures (4 and 9 weeks of treatment). Flash frozen leaves were individually homogenized 487

using a TissueLyser (Qiagen Retsch) and total RNA was isolated (from each leaf) using RNeasy Plant Mini 488

Kit (Qiagen) following the manufacturer’s instructions. The purity and integrity of total RNA extracts was 489

determined using a NanoDrop 8000 UV-Vis Spectrophotometer (Thermo Scientific) and 2100 Bioanalyzer 490

(Agilent), respectively. For each time point, RNA extracts from five leaves sampled from five different 491

plants were pooled and sequenced as a single sample. In addition, replicates from single individual leaves 492

were sequenced for selected timepoints (see Table S1 and “Differential expression” below). Two time 493

points lacked expression values: W9 in B. distachyon (RNA integrity was insufficient for RNA sequencing) 494

and W0 in S. lagascae (insufficient supply of plant material). Samples were sent to the Norwegian 495

Sequencing Centre, where strand-specific cDNA libraries were prepared and sequenced (paired-end) on an 496

Illumina HiSeq 2000 system. The raw reads are available in the ArrayExpress database 497

(www.ebi.ac.uk/arrayexpress) under accession number E-MTAB-5300.

498

Transcriptome assembly and ortholog inference

499

Using Trimmomatic v0.32 (Bolger et al., 2014), all reads were trimmed to a length of 120 bp, Illumina 500

TruSeq adapters were removed from the raw reads, low quality bases were trimmed using a sliding window 501

of 40 bp and an average quality cut-off of 15, and reads below a minimum length of 36 bp were discarded.

502

Read quality was checked using fastqc v0.11.2. For each species, transcripts were assembled de novo with 503

Trinity v2.0.6 (Grabherr et al., 2011) (strand specific option, otherwise default parameters) using reads from 504

all samples. We assessed the completeness of the transcriptomes using Benchmarking of Universal Single- 505

Copy Orthologs (BUSCO) v3.0.2 (Waterhouse et al., 2018) with the provided Embryophyta database, based 506

on OrthoDB version 9 (Zdobnov et al., 2017) (Table S6). Coding sequences (CDS) were identified using 507

TransDecoder rel16JAN2014 (Haas et al., 2013). Where Trinity reported multiple isoforms, only the 508

longest CDS was retained. Ortholog groups (OGs) were constructed from the five de novo transcriptomes 509

and public reference transcriptomes of H. vulgare (barley_HighConf_genes_MIPS_23Mar12), B.

510

distachyon (brachypodium v1.2), O. sativa (rap2), Z. mays (ZmB73_5a_WGS), S. bicolor (sorghum 1.4) 511

and L. perenne (GenBank TSA accession GAYX01000000) using OrthoMCL v2.0.9 (Li et al., 2003). All 512

reference sequences except L. perenne were downloaded from http://pgsb.helmholtz- 513

muenchen.de/plant/plantsdb.jsp.

514

14

High confidence ortholog groups

515

To compare gene expression across Pooideae, we required ortholog groups to contain one gene from each 516

species that all descended from a single gene in the Pooideae ancestor. As the ortholog groups (OGs) 517

inferred using orthoMCL sometimes cluster more distantly related homologs as well as include both 518

paraphyletic and monophyletic paralogs, we further refined the OGs by phylogenetic analysis. Several 519

approaches to phylogenetic refinement have been proposed previously (see e.g. Yang and Smith, 2014).

520

Here we first aligned protein sequences within each OG using mafft v7.130 (Katoh and Standley, 2013) 521

and then converted these protein alignments to codon alignments using pal2nal v14 (Suyama et al., 2006).

522

Gene trees were then constructed from the codon alignments using Phangorn v1.99.14 (Schliep, 2011) 523

(maximum likelihood GTR+I+G). Trees with apparent duplication events before the most recent common 524

ancestor of the included species were split into several trees. This was accomplished by identifying in-group 525

(Pooideae) and out-group (Z. mays, S. bicolor and O. sativa) clades in each tree, and then splitting the trees 526

so that each resulting sub-tree contained a single out-group and a single in-group clade. Finally, we only 527

retained the trees where all species in the tree formed one clade each (i.e. only monophyletic paralogs), B.

528

distachyon and H. vulgare formed a clade and at least three of the five studied species were included. These 529

trees constituted the high confidence ortholog groups (HCOGs).

530

Differential expression

531

Reads were mapped to the de novo transcriptomes using bowtie v1.1.2 (Langmead et al., 2009), and read 532

counts were calculated with RSEM v1.2.9 (Li et al., 2011). We observed higher numbers of monophyletic 533

species-specific paralogs in the de novo assembled transcriptome than in the reference genomes of H.

534

vulgare and B. distachyon. Since the de novo assembly procedure thus seemed to overestimate the number 535

of paralogs (most likely because alleles or alternative transcript isoforms were assembled into separate 536

contigs), we chose to represent each species in each HCOG by a single read-count value equal to the sum 537

of the expression of all assembled paralogs (analogous to so called monophyly masking (Smith et al., 538

2011)). By additionally setting counts for missing orthologs to zero, we created a single cross species 539

expression matrix with HCOGs as rows and samples as columns (Table S3).

540

To identify conserved and diverged cold response across species, differentially expressed genes (DEGs) 541

were identified using DESeq2 v1.6.3 (Love et al., 2014) with a linear regression model that combined the 542

species factor and the time point factor (with time points W4/9 as a single level). Pooled samples provided 543

robust estimates of the mean expression in each time point. To also obtain robust estimates of the variance, 544

the model assumed common variance across all time points and species within each HCOG, thus taking 545

advantage of both biological replicates available for individual time points within species and the 546

replication provided by analyzing several species. For each species, we tested the expression difference 547

between D0 and D1 (short-term response) and the difference between W0 and W4/9 (long-term response) 548

(Fig. S1). B. distachyon lacked the W9 samples and long-term response was therefore based on W4 only.

549

S. lagascae lacked the W0 sample and long-term response was therefore calculated based on D0. Due to 550

the observed diurnal effect (Fig. S2A), this might have resulted in more unreliable estimates of the long- 551

term cold response in S. lagascae. Genes with a false discovery rate (FDR) adjusted p-value < 0.05 and a 552

fold change > 2 or < 0.5 were classified as differentially expressed.

553

Sample clustering

554

Sample clustering was based on read counts normalized using the variance-stabilizing transformation 555

(VST) implemented in DESeq2. The purpose of the VST-transformation is to create expression values 556

where the variability is not related to the mean. For large counts, these VST-values are equal to log2- 557

transformed values. HCOGs that lacked orthologs from any of the five species, or that contained orthologs 558

with low expression (VST < 3), were removed, resulting in 4981 HCOGs used for the clustering. To

559

15 highlight the effect of the cold treatment over the effect of expression level differences between species, 560

the expression values were normalized per gene and species: First, one expression value was obtained per 561

time point per gene by taking the mean of the replicates. Then, these expression values were centered by 562

subtracting the mean expression of all time point. Distances between all pairs of samples were calculated 563

as the sum of absolute expression difference between orthologs in the 4981 HCOGs (i.e. manhattan 564

distance). The sample clustering tree (Figure S2A) was generated using neighbor-joining (Saitou and Nei 565

1987).

566

Comparison with known cold responsive genes

567

Genes identified as short-term cold responsive in H. vulgare were acquired from supplementary table S10 568

in (Greenup et al., 2011). These genes were found to be responsive to cold in three independent experiments 569

with Plexdb accessions BB64 (Svensson et al., 2006), BB81 (no publication) and BB94 (Greenup et al., 570

2011). The probesets of the Affymetrix Barley1 GeneChip microarray used in these studies were blasted 571

(blastx) against all protein sequences in our OGs. Each probe was assigned to the OG with the best match 572

in the H. vulgare reference. If several probes were assigned to the same OG, only the probe with the best 573

hit was retained. Correspondingly, if a probe matched several paralogs within the same OG, only the best 574

match was retained. DESeq2 was used to identify short-term response DEGs for all transcripts in all OGs 575

(i.e. this analysis was not restricted to the HCOGs), and these were compared to DEGs from (Greenup et 576

al., 2011). The statistical significance of the overlap between our results and those reported in (Greenup et 577

al., 2011) was assessed for each species by counting the number of genes that had the same response (up- 578

or down-regulated DEGs) and comparing that to a null distribution. The null distribution was obtained from 579

equivalent counts obtained from 100 000 trials where genes were randomly selected from all expressed 580

genes (mean read count > 10) with an ortholog in H. vulgare.

581

Gene ontology enrichment tests

582

Gene Ontology (GO) annotations for B. distachyon were downloaded from Ensembl Plants Biomart and 583

assigned to the HCOGs. The TopGO v2.18.0 R package (Alexa and Rahnenfuhrer, 2010) was used to 584

calculate statistically significant enrichments (Fisher’s exact test, P < 0.05) of GO biological process 585

annotations in each set of branch specific DEGs using all annotated HCOGs as the background. The 586

enrichment analysis was restricted to a subset of GO terms relevant to plants (i.e. plant GO plant). Branch 587

specific DEGs were those genes that were exclusively differentially expressed in all species within a clade 588

in the phylogenetic tree.

589

Positive selection tests

590

To test for positive selection, we first identified a species tree topology. HCOGs with a single ortholog 591

from each of the five Pooideae species and O. sativa were used to infer gene trees. BEAST v1.7.5 592

(Drummond and Rambaut, 2007) was run with an HKY + Γ nucleotide substitution model using an 593

uncorrelated lognormal relaxed clock model. A Yule process (birth only) was used as prior for the tree and 594

monophyly of the Pooideae was constrained. MCMC analyses were run for 10 million generations and 595

parameters were sampled every 10.000 generation. For each gene tree analysis, the first 10 percent of the 596

estimated trees were discarded and the remaining trees were summarized to a maximum clade credibility 597

(MCC) tree using TreeAnnotator v1.7.5. The four most supported gene tree topologies are shown in Fig.

598

S4. The placement of the lineages Meliceae and Stipeae within the Pooideae species tree is uncertain. The 599

species tree topology used in the positive selection tests was equal to the most common topologies among 600

the 3914 MCC trees.

601

Each of the HCOG gene trees were tested for positive selection using the branch-site model in codeml, 602

which is part of PAML v4.7 (Yang, 2007). We only tested branches for positive selection in HCOGs

603