Circadian and environmental signal integration in a natural population of Arabidopsis

Significance A major question in circadian biology is of how the molecular processes associated with circadian clocks operate in a natural setting. Using a natural population of Arabidopsis plants as an experimental model, we investigated how circadian timing signals are combined with information about the ambient light and temperature conditions, to regulate an intracellular signal transduction pathway. We identified multiple locations of environmental input within the signaling pathway, mathematically inferred causality within the pathway, and revealed 24 h fluctuations in environmental sensitivity, under field conditions, that are reminiscent of the process of circadian gating. Our approaches are of universal value for investigating circadian-regulated processes under natural conditions.


Supporting Text Figures S1 to S22 Tables S1 to S3 SI References
Other supporting materials for this manuscript include the following: Datasets S1 conditions across the field site, and sampled leaves from the locations nominated as "sun" and "shade" sites on successive days.At "sun" locations, plants received direct sunlight during the day, and at "shade" locations plants received sunlight filtered by surrounding vegetation for most of the day with the sites identified by measurement of the ratio of red to far red light (SI Appendix, Fig. S4; R:FR calculated as the photon irradiance from 660 to 670nm divided by the photon irradiance from 725 to 735nm (9)).During March 2015, plants received more direct sunlight, whereas during September 2015 the light was scattered through sky overcast with clouds in combination with persistent rain.
The second sets of samples were obtained under natural conditions with manipulation of the temperature conditions around patches of plants during 13 -14 September 2016, which was close to the autumn equinox at the field site.In addition to control plants that were not manipulated (SI Appendix, Fig. S14A), we applied two temperature treatments.These were (1) a continuous temperature increase (SI Appendix, Fig. S14B), whereby plants were covered with clear plastic horticultural domes to block air currents and trap warm air; (2) a continuous temperature reduction, using a custom device that passed air through a duct within a heatexchanging ice-filled polystyrene box and expelled the chilled air into a clear horticultural dome covering the plants, with chilling augmented by small ice packs within the dome (SI Appendix, Fig. S14C).

Field sampling for transcript analysis
Across all experimental conditions, the same sampling and RNA isolation procedures were used.
At 2-hour intervals, a fully expanded rosette leaf was excised with dissecting scissors from 6 replicate plants for each condition.The time-courses using naturally occurring sun and shade conditions each comprised 13 sampling timepoints over a total of 26 hours (from 14:00 on the first day to 16:00 on the second day).The time-courses involving artificial temperature manipulations comprised 15 sampling timepoints over a total of 30 hours (from 17:00 on the first day to 21:00 on the second day).The same replicate plants were sampled repeatedly through each time-series, but different plant patches were used in different sampling seasons.Sampled leaves were placed immediately into individual microtubes containing at least 400 µL RNAlater Stabilization Solution (Thermo Fisher Scientific, Waltham, MA, USA).Scissors and forceps were cleaned with 70% (w/v) ethanol between samples.After sampling, tubes were placed temporarily on dry ice for up to 2 hours, at -40 °C for 3 days in a portable freezer during transfer to the laboratory, and then at -80 °C until RNA isolation.During hours of darkness, sampling occurred using green-filtered head torches.
We wished to ensure that the abundance of transcripts could be compared between each sampling season.We normalized all transcript measurements to the transcript levels in one sample.Therefore, we obtained this reference sample for normalization of all RT-qPCR experiments in the study by pooling RNA from 10 leaves sampled at noon on 26 March 2015, from 10 healthy plants across the field site that were each separated by at least 1 metre.This provided a reference cDNA sample against which all RT-qPCR analyses from all sampling seasons were normalized within the DDCt method (10), to allow comparability between all datasets.The reference sample was collected at midday because all transcripts under investigation were expressed to some extent at that time point.In all experiments, dawn and dusk were defined as the astronomical (solar) time of sunrise and sunset.

RNA isolation and RT-qPCR
Frozen samples containing RNAlater were defrosted in a cold room for 4 hours, the RNAlater was removed, and leaf tissue was transferred to new dry tubes and frozen in liquid nitrogen.Frozen tissue was ground with a TissueLyzer (Qiagen, Hilden, Germany) and total RNA was isolated from the powdered plant material using Macherey-Nagel Nucleospin II RNA extraction kits (Thermo Fisher Scientific).RNA concentrations were determined using a Nanodrop spectrophotometer (Thermo Fisher Scientific).cDNA was synthesized using a High Capacity cDNA Reverse Transcription Kit (Thermo Fisher Scientific) and random primers supplemented with RNAase inhibitor (Thermo Fisher Scientific), as described previously (11,12).1:500 cDNA dilutions were analyzed using Brilliant III Ultra-Fast SYBR Green QPCR Master Mix (Agilent Technologies, Santa Clara, CA, USA), required primer pairs (Table S2), and Agilent Mx3005P qPCR instrument.Primers were designed using the PrimerQuest™ Tool (Integrated DNA Technologies, Coralville, IA, USA).Results were normalized using the DDCt method to AhgACTIN2 (11,12).AhgACTIN2 is encoded in A. halleri by locus g21632 (5) and has 97.8% coding sequence identity with A. thaliana ACTIN2 (At3g18780).We selected AhgACTIN2 as a reference transcript, because it has been used previously as a reference transcript for experiments involving A. halleri at this field site (2), AhgACTIN2 does not fluctuate across the seasonal cycle (1), and in our experiments, the RT-qPCR Ct for AhgACTIN2 did not oscillate across the diel cycle (AhgACTIN2 time-series arrhythmic after analysis with JTK_CYCLE test for rhythmicity (13); SI Appendix, Fig. S22).In time-series of observed data and model predictions, JTK_CYCLE was used to test for rhythmicity and estimate the phase, with phase expressed as the circular mean (calculated using the circular R package).Note that JTK_CYCLE reports the phase as an integer value.

Environmental monitoring
The temperature and irradiance were measured beside the plants during sampling.The temperature at each location, for each temperature manipulation, was monitored with EL-USB-2 data loggers (Lascar Electronics, Whiteparish, UK) at 5-minute intervals.Temperature loggers were wrapped in aluminium foil to prevent surface heating by solar radiation.Irradiance was measured using a CC-3-UV-S cosine corrector connected to a USB2000+ spectrometer with a QP400-2-UV-VIS fibre optic cable (Ocean Optics, Dunedin, FL, USA).Ambient light spectra (200 nm to 900 nm) were collected every 5 minutes over the 14 hours of light during each day of sampling using OceanView software (Ocean Optics) on a laptop PC, controlled by a custom script.The spectrometer and computer were powered using portable lithium battery packs (Powertraveller, Hampshire, UK).

Smooth trend model
The smooth trend model (STM) allows inference of a trend within time-series data that contains both sampling noise and biological stochasticity, and also allows a level of statistical confidence to be applied to that trend.The STM to assess the difference in transcript abundance between March and September under sun and shade conditions (Fig. 2) was defined by the equations: where  ! and  % are the smooth trend components in March and September in 2015, respectively,  is the time-varying difference between the two seasons,  ! and  % are the observed transcript abundance in the two seasons, and  % is the variance. = (1, 2, ⋯ , 13) is the time point at two-hour intervals.The same STM was used to analyze the difference in transcript abundance between sun and shade conditions in March and September (SI Appendix, Fig. S9).
For the models of the three (ambient, warm and cool) conditions in the temperature manipulation experiment (Fig. 4), additional ,  and  were considered: The parameters of the models were estimated by Bayesian inference using the Markov Chain Monte Carlo (MCMC) approach.We note that unlike classical hypothesis testing methods, multiple comparisons do not raise a problem in a Bayesian multilevel modelling (14).The statistical models were written in the Stan language and the programs were compiled using CmdStan (v2.24).To operate CmdStan, the cmdstanr package (v0.4.0) of R was used.After 1,000 warm-up steps, 1,000 MCMC samples were obtained by thinning out 3,000 MCMC samples for each of four parallel chains.Thus, 4,000 MCMC samples were obtained in total.We confirmed the convergence of MCMC sampling (SI Appendix, Figs.S5-S7).For  !,! ,  !,% ,  !, and  %,! , a flat prior over the entirety of real numbers (from negative infinity to infinity) was used to minimize prior influence on the posterior distributions.For , a flat prior over the range from zero to infinity was used, reflecting the underlying assumption that standard deviations cannot logically take negative values.A noninformative uniform prior is the first choice for standard deviations in hierarchical models because it does not constrain posterior inference (15), and it is also used in state-space models (16,17).We assumed that  ( (observation error) is common among  !,  % , and  ) , because it is mainly derived from technical errors during RT-qPCR which is assumed to be common among different samples.We confirmed the convergence of MCMC sampling for  ( (SI Appendix, Figs.S5-S7).A cauchy distribution was used to represent the time-varying difference between two time series () because this distribution corrects better for the influence of outliers, due to its relatively long tails and its efficiency for detecting change points in time series data (18,19).

Residual analysis for smooth trend model
Residuals were calculated as the difference between the mean observed values and predictions of the smooth trend model (STM), to assess the suitability of the model (SI Appendix, Fig. S8).
Quantile-Quantile plots of the residuals were drawn to compare its distribution to the normal distribution, using the qqplotr package (v0.0.5) of R.

Local level model with exogenous variables
The local level model with exogenous variables (LLMX) to analyze a diel trend and the effect of environmental variables on transcript abundance (Fig. 3) was defined by the equations: where  is the autoregressive trend component that is common among conditions,  is the regression coefficient,  is the true state of transcript abundance,  is the observed transcript abundance, and  % is the variance.The subscripts, , , , , ,  and ℎ The parameters of the models were estimated by Bayesian inference using the Markov Chain Monte Carlo (MCMC) approach (Table S2, S3).The statistical models were written in the Stan language and the programs were compiled using CmdStan (v2.24).To operate CmdStan, the cmdstanr package (v0.4.0) of R was used.After 3,000 warm-up steps, 1,000 MCMC samples were obtained for each of the four parallel chains, and thus 4,000 MCMC samples were obtained in total.We confirmed the convergence of MCMC sampling (SI Appendix, Figs.S10 and S12).
For  and  !, a flat prior over the entire real numbers (from negative infinity to infinity) was used to minimize prior influence on the posterior distributions.For , a flat prior over the range from zero to infinity was used, reflecting the underlying assumption that standard deviations cannot logically take negative values.A noninformative uniform prior is the first choice for standard deviations in hierarchical models because it does not constrain posterior inference (15), and it is also used in state-space models (16,17).We assumed that  ( (observation error) is common among  *+,-./,  *+,-6+40 ,  -02-./, and  -02-6+40 , because it is mainly derived from technical errors during RT-qPCR which is assumed to be common among different samples.We confirmed the convergence of MCMC sampling for  ( (SI Appendix, Figs.S10 and S12).

Embedding dimension
In the empirical dynamic modelling (EDM), time series is embedded into time-lagged series, a procedure known as state space reconstruction (20,21).The embedding dimension E is the dimension (i.e. the number of time-lagged series) used to reconstruct the state space.Prior to convergent cross mapping (CCM), the optimal E value was determined for each of AhgCCA1, AhgSIG5 and AhgpsbD BLRP in temperature manipulation experiments in September 2016, by univariate simplex projection (22), using the simplex function of the rEDM package (v0.7.5) of R.
We determined the optimal E values showing the maximum forecast skill ρ (SI Appendix, Fig. S17A).The optimal E value of 8 for AhgCCA1 was relatively large, considering those of AhgSIG5 and AhgpsbD BLRP, and data points available in this analysis.We adopted E value of 2 as the optimal value for AhgCCA1, because ρ is higher than 0.9 for E values from 2 to 8 (SI Appendix, Fig. S17A) and the model prediction well fits the observed data at E value of 2 (SI Appendix, Fig. S17B).The optimal E values determined were 2 for AhgCCA1, 4 for AhgSIG5, and 2 for AhgpsbD BLRP.The model predictions well fit the observed data for these genes at the determined E values (SI Appendix, Fig. S17B-D).

Convergent cross mapping
When a causal relationship exists from a variable X to another variable Y, the information of X can be found in the time series of Y, meaning that the prediction of X is possible using the information of Y.In convergent cross mapping (CCM) based on simplex projection, one predicts the nearest neighbors of X at time t in the reconstructed lagged trajectory, using their timecorresponding points of Y (23).When this prediction (called cross mapping) is successful, the causality from X to Y is assumed.The cross map skill (i.e.prediction skill) is evaluated by Pearson's correlation coefficient (ρ) between predicted and observed values.We performed CCM between AhgCCA1 (X) and AhgSIG5 (Y), and between AhgSIG5 (X) and AhgpsbD BLRP (Y) for the temperature manipulation experiments (September 2016), using the ccm function in the rEDM package (v0.7.5) of R. We considered time delay in the interactions by changing the parameter tp from −6 to 2 (negative, past; positive, future) in the ccm function.When causality exists, optimal predictability is expected to occur for tp ≤ 0, i.e, prediction of past values of X from Y (24).We used a technique known as multispatial CCM (25), to utilize the relatively small number of time points of 13/condition in 2015 and 15/condition in 2016.For the 2015 data, we used two light conditions (sun and shade) for each of March and September, giving a total of 26 time points applied to CCM.For the 2016 data, we used three temperature conditions (ambient, warm and cool), giving a total of 45 time points applied to CCM.To test the significance of cross map skill, we produced 1,000 diel surrogate time series for X into which a similar level of oscillation was incorporated, but the sequence of data deviation from the oscillation was randomized.We used three criteria for significant causality.First, optimal predictability occurs for tp ≤ 0, second, cross map skill is above the 95 % interval of diel surrogates at the optimal time lag (tp), third, cross map skill is improved according to the increase in a library size (number of time points used to reconstruct a state space), known as convergence (23).

Application of the LLMX model to the RNA-seq data
We       In each panel, the upper graphs show the predicted relative transcript abundance for sun (orange;   , equation 1 in Materials and methods) and shade (light grey;   , equation 3) conditions with the mean of observed values (dots), and the lower graphs represent the differences in transcript abundance between sun and shade conditions (, equation 2).The solid line and the shaded region are the median and the 95% confidence interval of the posterior distribution.When the 95% confidence interval of the difference between sun and shade conditions does not contain zero, the difference is considered significant and is indicated by asterisks.Panels E-J include estimation of rhythmicity and peak time relative to solar dawn of underlying data, using JTK_CYCLE.Vertical grey lines on time-series plots indicate the times of sunrise and sunset.STM analysis used data from 6 replicate plants per condition.for each time lag with a significant cross map skill.In B-F, the solid line represents the cross map skill between pathway components, and shaded area represents the 95% interval of the cross map skill using 1,000 diel surrogate time series as the explanatory variable that reflects the same degree of oscillation but with the sequence of variation randomized (i.e., null model).Cross map skill (ρ) provides a measure of the potential causality strength between the two variables.(D-F) Test of convergence, i.e., an improvement in cross map skill according to increase in a library size (number of time points used to reconstruct a state space), for each time lag with a high cross map skill.In B-F, the solid line represents the cross map skill between pathway components, and shaded area represents the 95% interval of the cross map skill using 1,000 diel surrogate time series as the explanatory variable that reflects the same degree of oscillation but with the sequence of variation randomized (i.e., null model).
represent temperature, irradiance, an upstream gene, March, September, sun condition and shade condition, respectively.The mean transcript abundance of the upstream gene (i.e., AhgCCA1 in the AhgSIG5 model and AhgSIG5 in the AhgpsbD BLRP model) was used as one of the inputs. = (1, 2, ⋯ , 13) is the time point at two-hour intervals.
used reanalyzed RNA-seq count data fromNishio et al. (2020) (26), with sequences originally obtained from four sets of 48-hour samples taken at 2-hour intervals in 2013(Nagano et al., 2019 Nat.Plants)(1).As environmental variables, we used ambient temperature at 2-hour intervals and hourly accumulated solar radiation recorded at the meteorological station in Osaka, the closest station to our field site with both environmental variables available.The same LLMX model (equations 9-17) was applied to each signalling pathway.For the modelling of AhgCOR15A (SI Appendix, Fig.S21), an additional term of a gene was added to equations 10-13.The parameters of the models were estimated by Bayesian inference as described in the previous section.

Fig. S1 .
Fig. S1.Relationship between AtCCA1 and AtSIG5 transcript abundance in A. thaliana under controlled conditions.(A-D) Relationship between AtCCA1 and AtSIG5 transcript abundance under conditions of constant light, from the transcriptome studies of (A) (27) (B) (28), (C) (29), (D) (30).(E, F) Relationship between AtCCA1 and AtSIG5 transcript abundance under light/dark cycles with (E) long and (F) short photoperiods, from the transcriptome study of (27, 31).Blue lines indicate a regression line.Pearson's correlation coefficient (R) with the corresponding pvalues, testing the likelihood of a chance correlation, are shown for each plot.

Fig. S2 .
Fig. S2.Data underlying models produced in this study, here comparing signalling pathway dynamics between two seasons, under different light conditions in a natural population of A. halleri.(A-D) Diel fluctuations in (A, B) ambient temperature and (C, D) total irradiance detected (200-900 nm), at 2-hour intervals (thinned out from original data measured at 5-minute intervals, for the purpose of aligning intervals with the transcript data) during sampling period in March and September 2015.(E-J) Transcript abundance of (E, F) AhgCCA1, (G, H) AhgSIG5 and (I, J) AhgpsbD BLRP.Vertical grey lines on time-series plots indicate the times of sunrise and sunset during March (solid line) and September (dashed).Panels E-J include estimation of rhythmicity and peak time relative to solar dawn of underlying data, using JTK_CYCLE.Data are mean ± s.e.m; n = 6 replicate plants.

Fig. S3 .
Fig. S3.Data underlying models produced in this study, here comparing signalling pathway dynamics between two different light conditions, during two sampling seasons in a natural population of A. halleri.(A-D) Diel fluctuations in (A, B) ambient temperature and (C, D) total irradiance detected (200-900 nm), at 2-hour intervals (thinned out from original data measured at 5-minute intervals, for the purpose of aligning intervals with the transcript data) during sampling period in March and September 2015.(E-J) Transcript abundance of (E, F) AhgCCA1, (G, H) AhgSIG5 and (I, J) AhgpsbD BLRP.Panels E-J include estimation of rhythmicity and peak time relative to solar dawn of underlying data, using JTK_CYCLE.Vertical grey lines on time-series plots indicate the times of sunrise and sunset.Data are mean ± s.e.m; n = 6 replicate plants.

Fig. S4 .
Fig. S4.The ratio of red to far-red light in a natural population of A. halleri, during March and September sampling seasons.(A, B) Comparison of the ratio of red to far-red light received by plants under the sun-and shade conditions during (A) March 2015 and (B) September 2015 sampling seasons.The R:FR varied during the photoperiod during both sampling seasons, and the effect of shade on R:FR was ameliorated by heavy cloud cover.Vertical grey lines on graphs indicate the times of sunrise and sunset.

Fig. S6 .
Fig. S6.Convergence of Markov Chain Monte Carlo (MCMC) sampling for smooth trend model (STM) for AhgSIG5 at the sun sampling site.(A) Density plots (left) and trace plots (right) of 1,000 MCMC samples/chain generated from posterior distributions of  (a smooth trend component),

Fig. S10 .
Fig. S10.Convergence of Markov Chain Monte Carlo (MCMC) sampling for the local level model with exogenous variables (LLMX) for AhgSIG5.(A) Density plots (left) and trace plots (right) of

Fig. S14 .
Fig. S14.Moderate temperature manipulations applied to patches of A. halleri plants, in the field, using custom-designed equipment.(A) Representative appearance of plant patches under naturally fluctuating conditions.(B) Plants covered with a plastic dome to raise the temperature.(C) Plants covered with a plastic dome undergoing temperature reduction with a custom chilling device.In this device, cool air is introduced to enclosed plant patches after being driven slowly through a heat exchanger, positioned within an expanded polystyrene box filled with ice.(D) Temperature changes in each condition at 2-hour intervals (thinned out from original data measured at 5-minute intervals, for the purpose of aligning intervals with the transcript data) during sampling period in September 2016.(E-G) Fluctuations in AhgCCA1, (F) AhgSIG5 and (G) AhgpsbD BLRP transcript abundance under ambient conditions and following temperature manipulation of plant patches.Panels E-G include estimation of rhythmicity and peak time relative to solar dawn of underlying data, using JTK_CYCLE.Grey shaded boxes on graphs indicate the period between sunset and sunrise.Data are mean ± s.e.m; n = 6 replicate plants.

Fig. S18 .
Fig. S18.Evaluation of causal relationships between components of a circadian signalling pathway in a natural population of A. halleri in September 2015.(A) The direction of causality (from X to Y) and that of prediction (cross mapping) in convergent cross mapping (CCM) (from Y to X) is opposite.(B-C) Estimation of causality between pairs of pathway components, across a range of time delays between the pathway component (time to prediction, tp) for data in September 2015.Cross map skill (ρ) provides a measure of the potential causality strength between the two variables.(D-F) Test of convergence, i.e., an improvement in cross map skill according to increase in a library size (number of time points used to reconstruct a state space), for each time lag with a significant cross map skill.In B-F, the solid line represents the cross map skill between pathway components, and shaded area represents the 95% interval of the cross map skill using 1,000 diel surrogate time series as the explanatory variable that reflects the same degree of oscillation but with the sequence of variation randomized (i.e., null model).

Fig. S19 .
Fig. S19.Evaluation of causal relationships between components of a circadian signalling pathway in a natural population of A. halleri in March 2015.(A) The direction of causality (from X to Y) and that of prediction (cross mapping) in convergent cross mapping (CCM) (from Y to X) is opposite.(B-C) Estimation of causality between pairs of pathway components, across a range of time delays between the pathway component (time to prediction, tp) for data in March 2015.Cross map skill (ρ) provides a measure of the potential causality strength between the two variables.(D-F) Test of convergence, i.e., an improvement in cross map skill according to increase in a library size (number of time points used to reconstruct a state space), for each time lag with a high cross map skill.In B-F, the solid line represents the cross map skill between pathway components, and shaded area represents the 95% interval of the cross map skill using 1,000 diel surrogate time series as the explanatory variable that reflects the same degree of oscillation but with the sequence of variation randomized (i.e., null model).

Fig. S20 .
Fig. S20.Prediction of AhgSIG5 signalling pathway dynamics in four seasons with different day length in a natural population of A. halleri.(A, B) Diel fluctuations in (A) ambient temperature and (B) solar radiation, at 2-hour intervals in the four sets of 48-h sampling period in 2013.(C) AhgCCA1 transcript abundance averaged over four biological replicates during the sampling period.(D) Bayesian estimation of the local level model with exogenous variables (LLMX) for AhgSIG5 transcript dynamics during the sampling period.Modelled transcript dynamics (lines and shaded area representing the median and the 95% confidence interval of the posterior distribution, respectively) are superimposed upon observed mean transcript abundance (circles).(E) Bayesian estimation of trend component of AhgSIG5 transcript (line and shaded area representing the median and the 95% confidence interval, respectively).Vertical dashed lines indicate the times of sunrise and sunset in each season.(F) Bayesian estimation of regression coefficient of environmental variables and a potential upstream regulator (AhgCCA1).Dots and error bars represent the median and the 95% confidence interval, respectively.Shaded boxes on A-D indicate the period between sunset and sunrise.LLMX analysis used data from 4 replicate plants per season.Analysis used RNA-seq data from (1).

Fig. S21 .
Fig. S21.Prediction of AhgCOR15A signalling pathway dynamics in four seasons with different day length in a natural population of A. halleri.(A, B) Diel fluctuations in (A) ambient temperature and (B) solar radiation, at 2-hour intervals in the four sets of 48-h sampling period in 2013.(C, D) (C) AhgCCA1 and (D) AhgCBF1 transcript abundance averaged over four biological replicates during the sampling period.(E) Bayesian estimation of the local level model with exogenous variables (LLMX) for AhgCOR15A transcript dynamics during the sampling period.Modelled transcript dynamics (lines and shaded area representing the median and the 95% confidence interval of the posterior distribution, respectively) are superimposed upon observed mean transcript abundance (circles).(F) Bayesian estimation of trend component of AhgCOR15A transcript (line and shaded area representing the median and the 95% confidence interval, respectively).Vertical dashed lines indicate the times of sunrise and sunset in each season.(G) Bayesian estimation of regression coefficient of environmental variables and potential upstream regulators (AhgCCA1 and AhgCBF1).Dots and error bars represent the median and the 95% confidence interval, respectively.Shaded boxes on A-E indicate the period between sunset and

Fig. S22 .
Fig. S22.AhgACT2 reference transcript was not rhythmic over diel timecourses in our experiments.Ct for AhgACT2 from RT-qPCR analysis of (A) sun and (B) shade datasets collected during September 2015.Data are mean Ct ± s.e.m; n = 6.p-values provided for statistical test for significant rhythmicity, using the JTK_CYCLE algorithm(13).

Table S1 .
Primers used for RT-qPCR analysis of transcript abundance in Arabidopsis halleri subsp.gemmifera.

Table S2 .
Estimated parameter values in the local level model with exogenous variables (LLMX) for AhgSIG5.The 2.5 %, 50.0 % and 97.5 % points of 4,000 MCMC samples obtained from posterior distributions are shown.

Table S3 .
Estimated parameter values in the local level model with exogenous variables (LLMX) for AhgpsbD BLRP.The 2.5 %, 50.0 % and 97.5 % points of 4,000 MCMC samples obtained from posterior distributions are shown.