Differential effects of day-night cues and the circadian clock on the barley transcriptome

The circadian clock is a complex transcriptional network that regulates gene expression in anticipation of the day-night cycle and controls agronomic traits in plants. However, in crops, information on the effects of the internal clock and day-night cues on the transcriptome is limited. We analysed the diel and circadian leaf transcriptomes in the barley cultivar Bowman and derived introgression lines carrying mutations in EARLY FLOWERING 3 (ELF3), LUX1, and EARLY MATURITY 7 (EAM7). Mutations in ELF3 and LUX1 abolished circadian transcriptome oscillations under constant conditions, whereas eam7 maintained oscillations of ≈30% of the circadian transcriptome. However, day-night cues fully restored transcript oscillations in all three mutants and thus compensated for a disrupted oscillator in the arrhythmic barley clock mutants elf3 and lux1. Nevertheless, elf3 but not lux1 affected the phase of the diel oscillating transcriptome and thus the integration of external cues into the clock. Using dynamical modelling, we predicted a structure of the barley circadian oscillator and interactions of its individual components with day-night cues. Our findings provide a valuable resource for exploring the function and output targets of the circadian clock and for further investigations into the diel and circadian control of the barley transcriptome.


116
Diel and circadian oscillations of the barley transcriptome 117 We analysed the diel and circadian global leaf transcriptome of the barley cultivar Bowman and 118 the derived introgression lines carrying mutations in HvELF3 (BW290), HvLUX1 (BW284) and 119 HvEAM7 (BW287). Plants were grown under cycles of 12h light and 12 h night (LD) and the 120 second leaf of replicate plants was harvested every four hours over 24h. Additional samples were 121 taken in a 2h interval at dusk in all genotypes and additionally at dawn in Bowman 122 (Supplemental Figure S1). Thereafter plants were transferred to constant light and temperature 123 conditions (LL) and leaf samples were taken every four hours for 36h starting from the first 124 subjective night. Individual libraries were single-end sequenced on a HiSeq 2500 with 10 125 Million reads per library and reads were mapped against a custom reference sequence consisting 126 of 68,739 transcripts (Digel et al., 2015). The nomenclature of the gene models used in this study 127 (Digel et al., 2015) were cross-referenced with the identifiers of the HORVU gene models 128 annotated on the barley pseudochromosomes (Mascher et al., 2017). Raw read counts normalized 129 to counts per million (CPM) were used for the downstream rhythmic analysis and modeling. We To investigate temporal expression patterns of the circadian-regulated transcripts under free-174 running conditions, we estimated the phase and the period of every circadian-regulated transcript 175 in the two genotypes that sustained free-running circadian rhythms, Bowman and eam7. In 176 Bowman, the distribution of the circadian transcriptome expression phase followed a bimodal 177 pattern with the highest number of transcripts peaking shortly before the transitions to subjective suggested that the free-running period under LL was extended in eam7 compared with Bowman.

189
Regulation of the transcriptome-wide phase in day/night cycles 190 Next, we investigated the transcriptome oscillations under the diel LD conditions. In all 191 genotypes, including those that were arrhythmic in LL, the mean of the period distribution was 192 consistent with the enforced 24-h diel cycle and ranged between 23.5 and 23.6 h (Supplemental 193 Figure S2). The phase was bimodally distributed over the day/night cycle in Bowman so that for 194 the highest number of transcripts the peak of expression occurred before dawn and dusk and, the  The analysis of the clock mutants, however, suggested that the bimodal phase distribution under 201 LD is controlled by both the circadian clock and day/night cues. In Hvelf3 the phase was 202 bimodally distributed under diel cycles similar to Bowman, however the quantitative 203 9 characteristics of the phase distribution differed. Namely, in Hvelf3, the phase distribution 204 showed higher peaks at dawn and dusk and deeper troughs during the night or the day than in   231 We then sought to infer the regulatory relationships between components of the barley circadian 232 clock. To this end, we modeled a transcriptional network based on the RNAseq time-series data.

233
Our data suggested that HvELF3 and HvLUX1 are integral components of the barley oscillator as 234 they were necessary to sustain transcriptome oscillations under LL ( Figure 1). Therefore, we 235 hypothesized that modeling a transcriptional network around HvELF3 and HvLUX1 could 236 identify the regulatory relationships that shape the circadian clock in barley. We followed an 237 approach that searches the dynamic dependencies of HvELF3 and HvLUX1 expression on other conditions. Therefore, only the expression datasets from Bowman and eam7 could be used in 243 modeling, since their transcriptomes oscillated under LL conditions. In both Bowman and eam7, 244 the transcripts encoding HvELF3 had a very low signal-to-noise ratio due to low rhythmicity 245 under LL and could not be used for modeling. We therefore rooted the network around HvLUX1, 246 which displayed robust oscillatory dynamics (Supplemental Figure S4). To reduce the 247 identification of erroneous interactions we filtered all circadian transcripts for those that were  We then investigated the consistency between the models obtained for Bowman and eam7 using 265 the v-gap metric (Supplemental Data 4, Supplemental Figure S5). This approach estimates 266 differences between models and allowed us to identify regulatory interactions that were 267 maintained or abolished in the eam7 mutant (Mombaerts et al., 2019). Following this approach,   (Faure et al. 2012, Campoli et al. 2013, Pankin et al. 2014.

496
BW290 carries an introgression of the early maturity 8 allele (eam8.k) which is characterized by 497 a base-pair mutation leading to a premature stop codon in HvELF3, which is orthologous to 498 ELF3 in Arabidopsis (Faure et al., 2012, Hicks et al. 2001  The quality of the sequencing data was verified using the FastQC software. The reads were 520 mappedagainst a custom barley reference transcriptome (Digel et al., 2015) and raw read counts 521 were obtained using the software implementing the full pipeline for RNA-seq analysis RobiNA 522 (version 1.2.3) with default settings (Lohse et al., 2012). Raw read counts were normalized to 523 counts per million (CPM) using the R/Bioconductor package "edgeR" and used for the 524 downstream rhythmic analysis and modeling. To cross-reference a nomenclature of the gene 525 models used in this study (Digel et al., 2015) with the identifiers of the HORVU gene models For the analysis of the day/night data, the sequence of samples was inverted to start with the 530 night followed by the day samples (12h, 14h, 16h, 20h, 0h, 2h, 4h, 8h), and this data set was WT and BW287 datasets, respectively. BW284 (Hvlux1) and BW290 (Hvelf3) datasets were not 558 considered in the following network analysis since these mutations led to the arrhythmic 559 transcriptomes. Finally, seven genes (Hv.21080, Hv.22191, Hv.23289, Hv.32914, Hv.33010, 560 Hv.6793, MLOC_7084.3) were manually discarded from both subsets list of candidates as they 561 were not DNA binding transcription factors but rather enzymes in a metabolic process, leaving We used first order models to represent the system dynamics between two genes at a time using 576 an LTI model with the following equation: With ( ) corresponding to the degradation rate of y and ( )that represents the influence of 580 another transcription factor through the synthesis rate of y. The model, therefore, evaluates 581 whether the rate of change of a particular gene y depends on another gene u. Estimating a model 582 means finding (a), (b) and (c) that produce a vector y(t) as close as possible to the real data. The 583 estimation of parameters was performed using the function 'pem' implemented in MATLAB that 584 minimizes the prediction error of the data. LL data were used for the estimation, as they 585 represent the autonomous behavior of the oscillator.

586
The goodness of fit of the model with the data was calculated as following:

Equation 2
588 589 Where is the data (output), ̅ is the average value of the data, and ^ is the estimated output.

590
MATLAB function compare was used to compute the fitness of the model. Each potential link 591 between two genes was validated if the associated model reproduced the dynamics involved with 592 a sufficient degree of precision, which corresponds to a fitness threshold estimated at 60% 593 (Supplemental Information).

595
To investigate the potential regulators of HvLUX1, a collection of independent 1 st order LTI 596 models was estimated separately between each of the transcript and HvLUX1 in the Bowman 597 background. In each case, the parameters were estimated so that they together provide the best 598 possible fit to the HvLUX1 time-course data. This step takes the following form:  Figure S4.

612
To further narrow down the predicted regulatory interactions, we estimated the consistency of 613 the candidate models using the filtered eam7 (BW287) dataset. For this purpose, we evaluated 1 st 614 order LTI models for each of the previously identified regulations and retained those with the 615 goodness of fit > 60% in the eam7 (BW287) experimental condition (Supplemental Figure S4).

616
To keep links with the highest confidence only, the dynamical consistency of the LTI models 617 based on these two independent datasets (Bowman and BW287) was evaluated using the nu (ν) Where ℎ was assumed to be binary (1 = light; 0 = dark). We fixed the light delay ℎ to 0h 649 to represent the effect of rapid light signaling on the transcripts, and computed delays ranging 650 from 0 to 8h, every 0.2h, for HvLHY. The delay that provided the best fit to the data was selected 651 independently for each transcript. Ultimately, models were validated if they succeeded in 652 capturing the regulatory dynamics involved with a goodness of fit > 60%. 653 We assessed the accuracy of our LTI-based network reconstruction algorithm on the circadian

Light-regulated
Phase Bowman LL (h) Figure 5: Relationship between external and internal cues to regulate the phase of the barley transcriptome a) Fractions of transcripts identified as clock-dominated, co-dominated by the clock and light and light-dominated by the Bode-analysis. b) Phase relationship between diel cycles (LD) and constant light (LL) for all transcripts oscillating in LD and LL and those dominated by the circadian clock, co-regulated by the circadian clock and light and light-dominated. c) Phase distribution of clock-dominated, co-regulated and light-dominated transcripts in diel cycles (LD).