Our experiment was designed to specifically assess the effect on rice global gene expression of climatic fluctuations, different types of field environments and the genetic background. We conducted two phases of field cultivation - one during the dry season (January–February) and one during the wet season (July–August) at the experimental rice station of the International Rice Research Institute (IRRI), Los Baños, Laguna, in the island of Luzon in the Philippines in 2013 (Figure 1A). The plants were cultivated in two adjacent fields typical of two different systems of rice cultivation. One field was cultivated following irrigated lowland practices: it was flooded with shallow water under constant irrigation and seedlings were transplanted to the field after being raised for 3 weeks in seedbeds (referred to as irrigated field). The second field was managed according to upland cultivation practices: it was rainfed, not irrigated and had been directly seeded (referred to as rainfed field). Each field was divided into two subfields constituting the biological replicates. We grew three different landraces of rice: (i) Azucena, an upland adapted landrace; (ii) Pandan Wangi, traditionally used for lowland cultivation only in the dry season; and (iii) Palawan, another upland adapted landrace that is used only during wet season cultivation.

Figure 1

Download asset Open asset

To avoid the major shift in gene expression patterns induced by the transition to the flowering stage (Sato et al., 2011), which would confound our detection of environmental effects, we sampled rice leaf tissue during 1 month of vegetative growth (15 sampling timepoints, 2 days apart). Sampling was carried out 4 hr after sunrise to minimize circadian-driven transcriptional variation. Each sample included six young leaves, each from a different plant, to minimize variation in individual plant microenvironment.

We measured global gene expression using RNA sequencing (Figure 1B). We excluded from our analysis genes for which we detected sequencing reads for less than 20 samples out of the 60 samples in each genotype per season subset of the data. The expression data for the 22,144 remaining genes was log-transformed and the biological replicates were averaged.

Our goal is to relate gene expression variation over time to variation in climatic conditions and plant developmental stage, and assess how these relationships are affected by season, field type and genetic background. We focused on trends in gene expression variation common to a high number of genes: after removing 1251 genes with a low coefficient of variation and 2962 genes with a low mean expression, we grouped the remaining genes into co-expressed gene clusters (Figure 1B). The number of clusters chosen was the highest that satisfied the constraint that no more than 5% of all the genes in the analysis belonged to “non-representative” small clusters, defined as containing less than 1% of all the genes in the analysis. We used the mean expression profile of all genes in each cluster as a representation of the variation in expression within that cluster.

We used a model selection approach to explain gene expression patterns by environmental and developmental variation. This approach relies on selecting a linear combination of environmental/developmental (ED) input parameters that both minimizes model mean squared error (MSE), quantifying the difference between the model and the expression data, and limits model complexity (i.e., avoiding over-fitting). A preliminary analysis showed that allowing for more than three parameters per equation over-fit the model more often than improving it, so we limited the number of parameters per linear equation to three. A typical ED equation had the following form: cluster mean = αED1 + βED2 + γED3, where ED1, ED2 and ED3 are ED parameters, and α, β and γ are linear regression coefficients. The ED parameters used in these models were measurements of current conditions at the time of sampling, recent changes in temperature, humidity, wind speed and solar radiation, temperature fluctuations, and short-term and long-term averages for all climatic conditions (Table 1). We included parameters for non-linear effect of short-term solar radiation on gene expression, because this type of effect has been observed on photosynthesis rate (Li et al., 2009). There were two parameters for field soil moisture, at 30 and 15 cm below ground (measured with tensiometers in the rainfed field and estimated to be constant at soil saturation value in the irrigated field), and a binary parameter for the field (irrigated or rainfed). A parameter was designed to represent developmental stage, using fixed values for the transplanting stage, end of tillering production and the heading time. Our set of ED parameters included those that correlated with each other, so we averaged nearly similar parameters (r > 0.98) and added to our model selection approach the constraint that two parameters with a Pearson correlation coefficient over 0.85 could not be selected in the same equation.

We designed our approach to take potential differences in climatic response between fields, genotypes and seasons into account. To assess whether disparities in transcriptional patterns could be explained using distinct ED equations, we considered different ED models for each cluster. The simplest model is a single equation for the whole cluster mean. A more complex model would combine two different equations, for example, in the case of a field-specific model: cluster mean = α_iED1 + β_iED2 + γ_iED3 in the irrigated field and cluster mean = α_rED4 + β_rED5 + γ_rED6 in the rainfed field.

If we were to take into account all possible differences between fields, genotypes and seasons, we would get models with as many as eight equations, which would be difficult to interpret. We therefore limited the maximum number of equations to four by applying the method to one season at a time, testing different equations between fields and genotypes; or considering only one genotype and examining the possibility of different equations between fields and seasons (this latter case is represented in Figure 1C). We chose between the models comprising one to four equations using the Bayesian Information Criterion (BIC), a statistical tool to limit over-fitting that includes a complexity penalty, calculated on the global model for all subsets.

We first used our model selection approach on the dry season data alone, as this was the season with the greatest phenotypic differences between the rainfed and irrigated fields. We wanted to determine the extent of genotype and field differences in gene expression variation during that season and how well the ED models can explain these differences. For a given cluster, we quantified the differences in variation of gene expression between the two genotypes by calculating the correlation between the two genotype-specific subsets of the cluster mean. The same method was used to evaluate differences in expression patterns between fields.

Only 19 out of the 56 clusters (4663 genes) had a correlation coefficient between genotypes below 0.8 (Figure 2A). Clusters with a low genotype correlation had a high model MSE, which showed that the ED models did not adequately explain these genotype differences. We found much more extensive differences in gene expression patterns between field environments (Figure 2B), with 48 clusters (15,103 genes) that had a correlation coefficient below 0.8. While all clusters with high correlation between field environments (r > 0.8) had ED models that fitted the cluster mean well (model MSE < 0.12), low MSE models were also selected for several clusters with strong dissimilarities between field environments. In some cases, ED models could thus explain strong field differences. We did not investigate the genotype effect further and instead used the genotypes as biological replicates that were averaged for the analysis of both seasons, incorporating all the expression data concatenated into 60 data points (Figure 1A).

Figure 2

Download asset Open asset

In this two-season analysis, in addition to identifying which ED parameters the expression of each gene cluster can be related to, we are assessing whether the field environment and the season affect the identified transcriptional response. As an example, the simple season-specific ED model selected for cluster 9 is presented in Figure 3A (Supplementary file 1A–C). The model consists of a single linear equation for both fields and both seasons. It combines a negative term for soil moisture at 15 cm depth, a positive term for 1-hr average of solar radiation (exponentially transformed) and a positive term for the change in temperature during the last 2 hr. In this model, the soil moisture input parameter (Figure 3B, Supplementary file 1B) allows for modeling gene expression differences between field environments while keeping a common model for both field environments, in particular the higher expression in the rainfed field environment in the dry season. The selection of a common model for the dry and wet seasons shows that the gene expression response to climatic factors of this cluster was largely independent of the season.

Figure 3

Download asset Open asset

We used genotype correlation within a cluster as a measure of replicability of ED effects: the higher the genotype correlation, the more gene expression response was driven by factors common to both genotypes (i.e., climatic and developmental factors). While the median genotype correlation of all genes in the analysis was 0.55, the median of all gene cluster means was 0.90, showing that averaging expression profiles over many genes remarkably reduces sources of non-replicability. We focused on 27 gene expression clusters with a genotype correlation greater than 0.9 (Figure 4, Supplementary file 1C). They encompass 11,371 (63%) of the 17,931 genes in the analysis, as they include most of the largest gene clusters. These clusters are also the ones with the best ED models (i.e., with a small error and low complexity, such as the one depicted in Figure 3, Supplementary file 1C). For these clusters, we observed extensive field effect on the climatic/developmental response, with only 6 out the 27 clusters (2308 genes) showing high correlation (r > 0.8) between fields in both seasons. As expected from the much lower precipitation levels during the dry season, which accentuated the difference in water availability between fields (Figure 1A), the dry season climatic conditions generated wider between-field differences in gene expression compared with the wet season (Figure 4, Supplementary file 1C). The season effect was slightly less prevalent, as 11 out the 27 clusters (4967 genes) were modeled with the same equation for both seasons (Supplementary file 1C), indicative of season-independent responses.

Figure 4

Download asset Open asset

We analyzed gene clusters according to their field response to understand how distinct modes of cultivation affect gene expression under the same climate. The expression of a gene cluster can be affected by the field environment in two ways: (i) distinct responses to climatic/developmental factors and/or (ii) a shift in expression level, representing different ways in which the effect of the field environment can be integrated with the climatic response and developmental program. Enrichment in specific functions or pathways within the different types of gene clusters can be indicative of the role of certain processes in the adaptation to distinct field environments.

We divided the 27-gene clusters into four groups based on whether they showed different expression responses to climatic/developmental factors in each of the two fields (correlation between the expression patterns of the two fields below 0.8) in one or two seasons (Figure 4). We also calculated the difference in mean expression between the fields for each season. To investigate the molecular processes that were most affected by climatic variation and field environments, we conducted a gene ontology (GO) term enrichment analysis for each cluster (Supplementary file 1E).

Correlations between ED parameters make it difficult to ascertain the causal factor of gene expression change from the parameter selected in a model. This is why, for the interpretation of ED models, we grouped those of the parameters that showed high correlations to each other (Figure 5A). The parameters (sometimes included in a group of highly correlated parameters) selected for the 27-gene clusters most representative of ED responses are shown in Figure 5B.

Figure 5

Download asset Open asset

The functional characterization of genes involved in the response to abiotic stresses is mostly conducted with single environmental perturbations under controlled conditions. Few genes have been functionally characterized for a role in the response to combined environmental stresses. We have found multiple clusters with patterns of expression showing additive effects of several environmental signals. Because changes in the expression of a gene in reaction to perturbations in an environmental factor are usually interpreted as this gene having a function in the physiological response to that factor, our results suggest the genes in these clusters are involved in the response to simultaneous environmental changes. We assessed this further by looking at whether genes shown to have a role in the response to a specific abiotic stimulus also responded transcriptionally to other environmental signals. We compiled lists of genes that had already been demonstrated to have a role in the response to light, water availability or temperature in rice as well as rice homologs of such genes in other species (mainly the model system A. thaliana). We looked at whether the expression patterns of these genes were consistent with their putative function. We focused our attention on the candidate genes that showed a replicable environmental response in our data (genotype correlation above 0.8) and belonged to one of the 22 gene clusters whose main ED model terms are environmental parameters (Supplementary file 1C), which we refer to as genes with a strong environmental regulation.

Light response genes

From several published reviews, we identified 273 genes involved in light response (Galvão and Fankhauser, 2015), in particular phytochrome signaling (Wang and Wang, 2015), phototropin signaling (Sullivan et al., 2009), retrograde signaling (Häusler et al., 2014; Vogel et al., 2014), UV signaling (Singh et al., 2014), brassinosteroid biosynthesis and signaling (Zhang et al., 2014) and the regulation of photosynthesis in response to light (Kaiser et al., 2015; Pfannschmidt and Yang, 2012; Rochaix et al., 2012). Out of the 45 genes with a strong environmental regulation (Figure 6A, Supplementary file 1F), 31 were in nine clusters modeled with positive terms for short-term temperature/solar radiation while 14 were in seven clusters modeled with negative terms for short-term temperature/solar radiation. Most of these genes were also affected by the field environment, especially in the dry season. However, one gene from cluster 15 was highly variable over time (it belonged to the 95^th percentile for variance among the 17,931 genes in the 53 clusters) but not between fields. It is an ortholog of A. thaliana ELIP1 and 2 genes, which are possibly involved in pigment accumulation in response to stress (Casazza et al., 2005; Rossini et al., 2006).

Figure 6

Download asset Open asset

A group of co-expressed genes included a zeta-carotene desaturase and a phytoene desaturase, both carotenoid biosynthesis enzymes (Fang et al., 2008), two FtsH proteases and an ortholog of the chloroplast sensor kinase (CSK) all of which have a potential role in photoprotection. OsCHS, OsC4H and orthologs of A. thaliana UVR2 and PDX1.3, four genes that could be involved in UV acclimation (Singh et al., 2014), were also co-expressed with the photoprotection genes. These results show that transcriptional regulation is an important component of light response, especially for photoprotection and UV acclimation, and integrates responses to other abiotic stimuli linked to the field environment.

Here we examine the correspondence of gene expression and climatic variation between our Philippine experiment and a previously published rice field transcriptomic study conducted in different (temperate) climatic conditions (Nagano et al., 2012). From the latter experiment conducted in Japan, we selected 52 daytime timepoints that spanned nearly 7 weeks of vegetative growth in irrigated fields for the Japonica cultivar Nipponbare. We centered (subtracted the mean) the data by time of day to eliminate circadian clock effects. We used the same climatic measurements as those found in our experiment, which were recorded every minute. We refer to this data as the partial Nagano dataset (PND). To determine to what extent our results would hold true under a different climate, we performed a new analysis of our data that excluded measurements from the rainfed field (using the irrigated field half of the two seasons analysis dataset), as the experiment in Japan was exclusively conducted in an irrigated field (Figure 1B). We used these results to test whether ED parameters selected in the models for our data could also explain expression variance of the same clusters of genes in the PND.

The models selected for the 60 gene cluster means calculated from our data (Supplementary file 1G) could be represented by a single equation for the irrigated field in both seasons (season-independent) or of distinct equations for the irrigated field in each season (season-specific). We reasoned that climatic response patterns that were not consistent between the two seasons of our experiment (season-specific models) were unlikely to be reproduced in the PND. Consequently, we only tested the transferability of models for those of our clusters with a season-independent model, when the model explained more than half of the variance of the cluster mean. For each of these 36 clusters, we calculated the mean expression profile in the PND of the genes in the cluster. We determined whether one or several parameters in the model of the cluster mean over our irrigated field data could model the cluster mean over the PND with the same coefficient signs, by comparing all possible models with combinations of these parameters with the BIC.

Non-null models were selected for 24 of the 36 PND cluster means (Figure 7). They explained between 13 and 86% of the expression variance over the 52 data points of the PND. The parameters that could explain expression variation in our experiment and the PND were predominantly developmental stage and parameters for short-term solar radiation and temperature, especially in the models that explained a large part of the expression variance. Only one long-term climatic parameter was conserved in several PND cluster mean models: the 15-day average of wind speed, which was selected in four models, including one that explained 52% of the expression variance. These results show that a large part of climatic effects on gene expression vary with the type of climate; in particular, long-term effects.

Figure 7

Download asset Open asset

As we found that several clusters enriched for genes involved in photosynthesis and development were among the most sensitive to the field context, we investigated whether this was also true for the seasonal and climatic context. We performed a GO term enrichment analysis on the clusters from the irrigated field analysis (Supplementary file 1I). We identified two clusters highly enriched for genes related to photosynthesis (p<10^–27 and 10^–26, respectively), which were both modeled with season-dependent models. This result indicates that the transcriptional response of groups of co-expressed genes involved in photosynthesis is affected by the seasonal context. Two clusters (13 and 25) were enriched for genes associated with the developmental process and cell cycle. Cluster 13 was modeled with a season-dependent model while cluster 25 could be modeled with the same equation for both seasons of our experiment but no parameter of this equation could be transferred to model expression in the PND. This shows that the transcriptional regulation of some genes involved in development is not only sensitive to the field context but also to the season and climate type context.

To compare our method with the method from Nagano et al. (2012) with respect to the detection of climatic effects, we conducted another analysis of the PND using our model selection approach, determining clusters and models independently of the results obtained on our data (Supplementary file 1H). About a half of the models determined with the method by Nagano et al. (2012) do not have an environmental term and those that do, include the effect of a single climatic factor per gene, with a term that can be non-linear with diurnal changes in sensitivity. This makes their models intrinsically incompatible for a systematic direct comparison with ours. We thus chose to focus on two sets of genes identified by each method as having a clear environmental response and no developmental effect, and then evaluate the results of the alternative method on these genes. First, we selected from the results of our analysis genes highly correlated with their respective cluster means (r > 0.8) and belonging to three clusters with a low error model (MSE < 0.25) that included neither a developmental stage term nor a climatic parameter strongly correlated with the developmental stage. Among these 390 genes, 156 (40%) had a Nagano et al. (2012) model that included an environmental term for the timepoints we considered (between 8:00 and 16:00). Second, we chose genes with Nagano et al. (2012) models, including no developmental or circadian terms and with a linear environmental term for either temperature or solar radiation (which have the most reproducible effects) that explained more than 50% of the expression variance. When looking at the results of our analysis for these 214 genes, we found that 155 (72%) of them were in a cluster with at least one environmental term and a model error lower than 0.35 (R² > 0.7). Our method was thus more efficient in matching the Nagano et al. models with a clear environmental term than the Nagano et al. (2012) method was in the reciprocal comparison. This indicates that our approach has a higher sensitivity to the effect of climatic factors on gene expression, in addition to its ability to detect the co-occurring effects of multiple environmental parameters.

We used a model selection approach to identify relationships between major variations in global gene expression and environmental conditions/developmental stage. Focusing our analysis on groups of genes showing consistent variation across two different genotypes in each of the two seasons, we determined that most of these representative expression patterns could be explained through the combined effects of several environmental parameters, related to distinct climatic/soil-related factors and/or on different time-scales.

Co-occuring abiotic stresses trigger complex responses that cannot be predicted from the effect of single stresses (Prasch and Sonnewald, 2015), especially at the level of transcriptional regulation. Our dry season experiment, which included a drought period, presented ideal conditions to understand how limited water availability is integrated with other climatic signals at the transcriptional level in the field in comparison to controlled conditions. One common finding of studies on the combined effects of drought and heat is that a high number of genes are differentially expressed only by the combination of both stresses (Johnson et al., 2014; Rizhsky et al., 2004, 2002b). However, these results might be due to the fact that the differential expression induced by individual stresses is under the threshold of significance, while the addition of both small effects passes the differential expression threshold. Part of the assumed specificity could therefore be an artifact of the analysis. Although we cannot make a direct comparison with our analysis, as we do not have conditions testing the effect of drought alone, we can contrast gene expression responses between irrigated and rainfed fields during the dry season to assess how drought impacts the response to other abiotic stimuli and test whether there is indeed a specific response to drought combined with short-term climatic signals. Excluding clusters from groups 3 and 4, the expression of which is impacted by non-drought related field effects in the wet season, the different effects induced by drought in the rainfed field for clusters modeled with terms for short-term averages of temperature/solar radiation in the irrigated field during the dry season were (i) show the same response to short-term temperature/solar radiation with no or little shift in expression (clusters 18, 21 and 24); (ii) show an additive effect of drought over the short-term climatic response (clusters 4, 9, 12 and 30); or (iii) not display the short-term climatic response under drought (clusters 1, 8, 10 and 11). Only cluster 44, containing 204 genes, responded to short-term temperature/solar radiation variations in the dry season, specifically in the rainfed field. Our results therefore show that the main behaviors under combined drought and short-term abiotic signals are an addition of both responses and response to either one of the stimuli (often integrated with a long-term climatic response), while responses specific to the combination of both stresses are rare.

Genes that have been previously characterized for their role in a specific type of abiotic stimulus – light, temperature or water availability – and with a replicable expression pattern in our experiment generally displayed a transcriptional response to that stimulus. In many cases, these genes also seemed to respond to other abiotic stimuli, demonstrating the extensive connections between environmental pathways occurring in complex natural conditions. Some genes did not, however, have a pattern of expression consistent with what was measured in previous laboratory studies. For example, OsCDPK7 and OsMIOX have been shown to be induced by drought (Duan et al., 2012; Wan et al., 2007). Nevertheless, in our data, there were no differences in mean expression between the irrigated and rainfed fields in the dry season for OsCDPK7, while OsMIOX was less expressed in the rainfed field than in the irrigated field during both seasons. These discrepancies may be explained by the fact that abiotic stresses are generally analyzed individually. Some genes induced by drought, when water availability is the only varying environmental parameter, might be mostly responsive to temperature during a combination of drought and heat stress.

Another possible explanation is the difference in time-scale of our field study compared with most laboratory-based experiments. The effect of an abiotic stress is often analyzed within hours of the stress treatment; in our field experiment, we examined climatic variables that changed constantly and gradually over a month. The time-scale of the study is especially critical in the case of drought stress, because in natural conditions, it usually occurs over the course of days or weeks. More generally, covering longer time frames than laboratory studies allowed us to observe environmental effects interact with plant developmental processes in an integrated way (Allahverdiyeva et al., 2015), as well as examine the effects of long-term seasonal climatic effects. Our results thus show that measuring gene expression under natural conditions and for long periods of time can help better assess the effect of multiple abiotic stimuli on gene expression.

In this study, we assessed the impact of different seasons, field environments and types of climate on weather-driven patterns of gene expression. Most of the seasonal differences were limited to the rainfed field and were related to the prolonged period, with very little precipitation occurring during the dry season. When analyzing only the irrigated field data, we found that 44 out 60 clusters (representing 72% of the genes in all clusters) were modeled with a common equation for both seasons (Supplementary file 1G). A fair level of reproducibility was also observed between different years of irrigated field culture for rice (Nagano et al., 2012).

In contrast, the type of field in which rice was cultivated had a major effect on the response to climatic conditions of many genes. A large part of the differences in climatic response between field environments that we detected were especially pronounced in the dry season and linked to the water status of the rainfed field. Some expression differences, however, seemed independent of this drought effect, as they were detectable in the wet season.

Finally, we used independent expression data to assess the effect of a different type of climate on transcriptional responses in irrigated conditions. We found that only a small subset of the gene expression patterns observed under the tropical climate of our experiment could be generalized to the temperate conditions of the Nagano et al. (2012) study. This result shows that we can expect to reproduce findings about gene expression despite seasonal variations, but we should be careful about generalizing results to other climates or soil conditions than the ones their study was set in.

When investigating the reproducibility of transcriptional responses, we found that gene expression patterns driven by short-term averages of solar radiation and temperature, and to a lesser extent long-term wind speed, were the most consistent responses across seasons, fields and climates. Short-term averages of solar radiation and temperature belong to the group of parameters that was used to model most of the 27 gene co-expression clusters, with the highest correlation between genotypes in our two-season analysis. In particular, the exponentially transformed short-term averages (15 min, 1 or 4 hr) of solar radiation (NL− ) that model stronger effects of variations in the lowest range of irradiance were selected to model at least one entire season or field subset for 12 of the 27 gene clusters of our analysis and could systematically be transferred from the models of our irrigated field data to explain the expression variation of the same clusters of genes observed under the Japanese temperate climate (Supplementary file 1G). In the conditions of the Nagano et al. (2012) experiment, short-term solar radiation variation was weakly correlated with variation in temperature and precipitation, showing that this effect is likely to be specific to irradiance level. The reason why variation in the lowest range of solar radiation values seem to have a stronger effect on gene expression (compared with higher irradiance levels) for some clusters needs to be investigated. One hypothesis is that responses to low light can rely mostly on changes in gene expression because they do not require as rapid a response as high light conditions (Pfannschmidt and Yang, 2012).

In contrast, we found that the effect of longer-term variation (10-, 15-day averages) in temperature and solar radiation were much less reproducible across seasons, fields and climates than short-term effects. One explanation is that, as these responses may be integrated over long periods of time, there is greater opportunity for them to be modulated by broad climatic and developmental factors. The only long-term environmental factor that could be transferred to model gene expression under the Japanese climate for several clusters was the 15-day average of wind speed. In our two-season analysis, it was selected to model across both seasons and fields cluster 32, which is enriched for thigmotropism genes. While there have already been studies aiming at identifying genes responsive to mechanical cues (Lee et al., 2005), they were limited to hour-scale laboratory experiments. Our results show that, in nature, the effect of wind on gene expression might be more prevalent on a week-scale, which is consistent with the long-term accommodation to repetitive wind loads observed in poplar (Martin et al., 2010). Further study of the genes in cluster 32 should contribute to the understanding of gene regulation in response to wind necessary to unravel the cellular processes involved in the time integration of mechanosensing (Moulia et al., 2011).

The analysis for GO annotation enrichment of the 27 clusters with the most reproducible expression patterns revealed functionally related groups of genes that are co-expressed in response to environmental and developmental signals. These results point to the biological processes most critical for plant physiological response to dynamic environments. The GO terms for biological processes with the most significant enrichment were nucleic acid metabolic process (and its parent term nucleobase-containing compound metabolic process) and photosynthesis.

Nucleic acid metabolic process related genes were highly enriched in clusters 4 and 6 (p<10^–18 and 10^–22, respectively), which were also both enriched for developmental process and cell cycle associated genes. The strong enrichment in nucleic acid metabolic process related genes highlights the importance of a tight transcriptional regulation of this molecular process, probably linked to the control of cell division, in the adaptation of the developmental program to environmental fluctuations (Rymen and Sugimoto, 2012). Both clusters had a negative correlation between field expression profiles in the dry season, while the field correlation was above 0.75 in the wet season, during which both cluster means were modeled with a positive term for the developmental stage parameter. This indicates an extensive impact of drought on the expression pattern of the genes in these clusters. As cluster 6 was modeled with several short-term environmental parameters, our results show that the transcriptional regulation of cell division-related aspects of growth integrates with the developmental program several time-scales of environmental stimuli.

Clusters 1, 5 and 13 were greatly enriched for photosynthesis related genes (p<10^–17, 10^–9 and 10^–19, respectively). These gene clusters were all down-regulated in response to drought, in accordance with previous results in A. thaliana (Chaves et al., 2009). Clusters 1, 5 and 13 were also modeled with short-term climatic parameters. These effects were detected mostly for the irrigated conditions, with positive terms for solar radiation in the case of clusters 1 and 5. In contrast, cluster 13 was modeled with a negative term for short-term temperature in both fields. These differences indicate that different components of the photosynthetic machinery might need to respond to different environmental stimuli. Although it has been thoroughly demonstrated that a large part of environmental effects on photosynthesis are regulated at the physiological and biochemical levels (Kaiser et al., 2015), the fine regulation of transcript abundance for these photosynthesis genes shown here confirms that the control of gene expression is also an important component of the modulation of photosynthetic activity (Pfannschmidt and Yang, 2012).

Focusing on our irrigated field data and the data from an independent field experiment (Nagano et al., 2012), we found that the expression patterns of groups of co-expressed genes involved in photosynthesis were affected by the seasonal context (independently of drought) and the environmental response of groups of genes associated with development was dependent on the climatic context. This suggests that the transcriptional regulation of both photosynthesis and development relies on complex mechanisms, resulting in the integration of numerous layers of environmental cues, short-term climatic fluctuations as well as steadier aspects of the plant surroundings.

Understanding the dynamics of gene expression is a major challenge in biology. This is particularly difficult in the context of complex environments found in nature, as it requires unraveling the concurrent effects of multiple, fluctuating environmental signals on transcriptional patterns. Nevertheless, the ability to establish the environmental response of whole transcriptomes can have wide applications, including controlling engineered gene circuits (Uhlendorf et al., 2012) and predicting gene patterns in untested conditions (Bonneau et al., 2007; Danziger et al., 2014; Nagano et al., 2012). Our results suggest that while field environments can result in complex responses, one can nevertheless identify co-expressed gene clusters the mean expression of which can be accurately modeled with climatic, field, seasonal and developmental factors. Further work can integrate such models with other approaches, including gene network inference (Bonneau et al., 2007), genotype-by-environment interactions (Marais et al., 2013) and phenotypic modeling (Aikawa et al., 2010; Satake et al., 2013), to provide a more comprehensive picture of plant responses in natural environments. This, in turn, can be used in the design of improved crops (Hammer et al., 2006; Mochida et al., 2015) and the prediction of the ecological effects of climate change (Stafford et al., 2013).

We used rice landraces for which seeds were available at the IRRI in Los Baños, Philippines, and that we knew are traditionally used for either lowland or upland cultivation. Azucena and Palawan are Filipino upland landraces while Pandan Wangi is a lowland Indonesian landrace. Azucena and Padan Wangi were grown in the dry season, while in the wet season we grew Azucena and Palawan.

The rice plants were grown at IRRI in a 25 × 90 m field divided into equal size rainfed and irrigated sub-fields. The dry season and wet season experiments took place during the months of January–February 2013 and July–August 2013, respectively. The two seasons differed for the amount and frequency of precipitation as well as for temperature, which is generally lower during the dry season. During our experiment, the mean daytime temperature was 26.8°C for the dry season and 27.6°C for the wet season. The total precipitation during the experiment period was 84.5 mm in the dry season and 427 mm in the wet season. Precipitation of more than 1 mm was observed for 6 out of the 29 days of the dry season experiment and 22 out of the 29 days of the wet season experiment (Supplementary file 2).

In the irrigated field, one 21 day-old seedling was transplanted per hill with a spacing of 20 × 30 cm. For the rainfed field, direct seeding was performed with four seeds per hole, with a spacing of 20 × 30 cm; seedlings were thinned out to one per hill after 2 weeks. The composition of the fertilizer used was 120:30:30 nitrogen/phosphorus/potassium (NPK); it was applied as recommended. Carbofuran insecticide/nematicide was applied at 1 and 15 days after sowing (DAS) and herbicide was applied at 1 DAS in the upland ecosystem. Manual weeding and general plant protection were performed as needed.

The first sampling took place 16 days after transplanting seedlings in the irrigated field in the dry season and 23 days after transplanting in the wet season. Each sample consisted of six young leaves (of approximately the same size throughout the experiment) from six individual plants. Each leaf was immediately frozen in liquid nitrogen upon collection. We tried to reduce as much as possible the effect of circadian variation on gene expression, first between sampling time-points, by always starting the collection exactly 4 hours after sunrise. Second, to avoid a shift in expression within a sampling time-point due to the delay between the first and last collected samples, we ensured fast sampling by marking beforehand each plant and leaf to be collected. Collection took, on average, 13 min (between 11 and 15 min, except for the first time-point of the dry season, which took 18 min).

Averages of temperature, relative humidity, wind speed, solar radiation, precipitation and atmospheric pressure at the field site were recorded every 15 min using Wireless Vantage Pro2 ISS with 24-hr fan aspirated radiation shield from Davis Instruments (Hayward, CA). Measurements started 15 days before the first sampling day. After the ninth sampling of the dry season and during all of the wet season, from 1 hr before sampling to the end of sampling, we recorded climatic averages every minute. There was a weather station in the irrigated and rainfed fields, one of which also included a solar radiation sensor. No significant difference in measurements was detected between the fields, so the data of the station with the solar radiation sensor was used. We used 60 cm long 2710ARL tensiometers from Soilmoisture Equipment Corp. (Goleta, CA) placed 30 and 15 cm deep in the soil at four different locations, two per replicate of the rainfed field.

Plant height and tiller number were measured every 6 days for the same plants all along each season until the end of sampling (four plants per replicate in the dry season, twelve plants per replicate in the wet season).

Frozen leaf tissue was ground manually with a mortar and pestle cooled in liquid nitrogen. Total RNA was extracted from about 200–400 µl of ground tissue using the RNeasy Plant Mini Kit (Qiagen, Venlo, The Netherlands), following the manufacturer’s protocol and eluting the RNA in 40 µl of water. RNA was treated with Baseline-ZERO DNase (Epicentre, Madison, Wisconsin) according to the manufacturer’s instructions then cleaned up with the Qiagen RNeasy Mini Kit and eluted in 32 µl of water. We assessed RNA quality using nanodrop and electrophoresis on an agarose gel. Total RNA, 4 µg, were depleted of ribosomal RNA using Epicentre Ribo-Zero Magnetic Kit for plant leaf tissue. We purified the depleted RNA with the Agencourt RNAClean XP kit (Beckman Coulter, Brea, CA). We constructed RNA libraries using the Epicentre ScriptSeq v2 RNA-Seq Library Preparation Kit, purifying the complementary DNA and libraries with the Agencourt AMPure XP System. We added barcode index using the Epicentre ScriptSeq Index PCR Primers. We quantified the libraries by Qubit (Life Technologies, Norwalk, Connecticut), with the DNA HS kit. Libraries quality and average fragment size was determined using the 2100 Bioanalyzer (Agilent, Santa Clara, CA) with high sensitivity DNA reagents and DNA chip. We quantified the libraries on the LightCycler 480 (Roche, Nutley, NJ) using the KAPA (Wilmongton, MA) Library Quantification Kit. Libraries were sequenced using HiSeq 2000 (Illumina, San Diego, CA) 51 bp paired-end sequencing, with either 6 or 8 libraries per lane. Each sample provided a mean of 58 million sequencing reads.

The reads were aligned to the O. sativa Nipponbare – release 7 of the MSU Rice Genome Annotation Project reference genome (Kawahara et al., 2013) – 373,245,519 base pairs of non-overlapping rice genome sequence from the 12 rice chromosomes. Also included are the sequences for chloroplast (134,525 bp), mitochondrion (490,520 bp), Syngenta pseudomolecule (592,136 bp), and the unanchored BAC pseudomolecule (633,585 bp). The annotation contains 56,143 genes (loci), of which 6457 had additional alternative splicing isoforms resulting in a total of 66,495 transcripts.

We used Tophat (Kim et al., 2013; Trapnell et al., 2009) version 2.0.6 to align the reads, discarding low-quality alignments (quality score below 1). To count the number of reads that uniquely mapped to genes we used HTSeq (Anders et al., 2015) version 0.6.1. We compensated for variable sequencing depth between samples using the median-of-ratios method of DESeq2 (Love et al., 2014) version 1.6.3, and further performed a variance stabilizing transformation provided by the same package. We used the normalized count data for downstream analysis. This data is available in the Gene Expression Omnibus database (http://www.ncbi.nlm.nih.gov/geo/) under the accession number GSE73609.

We conducted a multi dimensional scaling of the normalized expression data where the samples clustered by genotype and field and to a lesser extent, season. We used these results to detect potential sample mislabeling and identified one sample switch. We also found that the sample for the first replicate of the Pandan Wangi rainfed field in the dry season, sixth time-point, did not cluster with any of the genotype/field groups so we removed it from the analysis and replaced it with a duplication of the second replicate.

We considered each subset of the data that consisted of one genotype and one season and excluded from the normalized expression dataset genes for which we detected no read for more than 40 samples in any of these subsets, which reduced the dataset to 22,144 expressed genes. We transformed the obtained value with the following function: log2(x + 1), as to keep positive values of expression and averaged the biological replicates. We calculated the coefficient of variation of the log-transformed expression over the 240 data points for each expressed gene. We identified 1251 genes with a coefficient of variation below 0.01, which were considered as having a stable expression in our experiment and were removed from the analysis. A preliminary cluster analysis showed that genes with a very low expression level had a weak correlation with the center of their clusters, thus 2962 genes with a mean lower than 1 were not included in the clustering step. To remove the absolute differences in expression level between genotypes and seasons, which we did not intend to model, each genotype/season subset was centered.

For each of the analyses, the expression data for each gene was scaled over the whole profile. Our clustering method was the Partitioning Around Medoids algorithm from the “cluster” package version 1.15.2 (Maechler et al., 2015) in R (R development Core Team, 2011) using 1 − r, with r the Pearson correlation coefficient, as the distance between expression profiles. To choose the number of clusters (k) we first used several clustering indices but they gave inconsistent results and usually indicated a k lower than 20. To avoid running the risk of under-clustering, which would have resulted in averaging genes responding to very different ED factors and therefore losing transcriptional signal, we chose to constrain k only on the fact that most clusters should represent a major trend in transcriptional variation. We therefore chose k as the highest number of clusters for which no more than 5% of all the genes in the analysis belonged to “non-representative” small clusters, defined as containing less than 1% of all the genes in the analysis. The gene co-expression clusters were ordered according to the number of genes they contained so that cluster 1 was the biggest cluster.

Two sets of gene annotations were obtained from Gene Ontology (http://geneontology.org/page/about) from the November 2014 release. For the first one, the annotation database was queried via MSU locus identifiers; for the second one, the database was queried via Uniprot Ids, obtained via a mapping from MSU Id to OMA Id, and then OMA Id to Uniprot Id (mapping files available in the current release of OMA). A third set of annotations was obtained directly from OMA (http://omabrowser.org/oma/about/). All three annotation sets were then combined non-redundantly in order to produce the final annotation file for rice genes. The enrichment analysis was conducted using the GOstats package in R (Falcon and Gentleman, 2007).

To represent each cluster, we calculated the average scaled expression of all the genes in the cluster. This “cluster mean” was used as the output data in our modeling approach. Climatic, field and developmental parameters constituted the input data. We calculated averages of temperature, relative humidity, wind speed, solar radiation, precipitation and pressure for 15 min 1, 4, and 24 hr; and 3, 6, 10 and 15 days before sampling. Short-term changes in temperature, humidity, wind speed and solar radiation were determined by calculating the difference between the value of each measurement at the sampling time and 20 min (δ20 min), 1 hr (δ1 hr) and 2 hr (δ2 hr) earlier, using 5, 10, and 30 min averages, respectively. We evaluated temperature fluctuations by decomposing daily variation in temperature with the seasonal decomposition by loess (stl) function in R and calculating 1 hr (ε1 hr), 4 hr (ε4 hr) and 24 hr (ε24 hr) before sampling averages of the remainder of the decomposition. The value used for the soil moisture parameters was a mean of measurements from four tensiometers, two in each replicate of the rainfed field. Exponentially transformed values of solar radiation used the following equations:

We designed a parameter that would give an estimate of the developmental stage, rather than age, of the plant to be able to compare appropriately the plants in the rainfed and irrigated fields as they followed different developmental itineraries relative to their age in days. This parameter used three stages as fixed points: transplanting stage (given a value of 0), corresponding to the actual transplanting event for the irrigated field and determined for the rainfed plants as the stage where they were the same height as just transplanted irrigated plants, end of tillering production (40) and heading stage (100). Intermediary time-points were calculated linearly between these fixed points. Input parameters were centered per genotype and season and scaled over the whole dataset to match the expression data.

When input parameters were more highly correlated than the highest genotype correlation of all clusters (r = 0.98), the correlated parameters were averaged together, as we would not have enough precision in the expression data to discriminate between them.

The same method was used on each cluster of a given set to select a model. It is described here for the analysis of one genotype in two seasons Source code 1; the analysis of two genotypes in one season is identical and the analysis of the irrigated field with averaged genotype follows the same principle.

1. We used the high domensional inference (hdi) function of the hdi R package (Meinshausen et al., 2009; Meinshausen and Bühlmann, 2010) using the lasso algorithm for the parameter selection, with the “stability” method, B = 300, EV = 2, threshold = 0.65 and fraction = 0.85, to select a stable subset of all the input parameters using the whole 60 data points.

2. All possible combinations of one, two and three parameters from this subset were used to calculate linear regression models, the MSEs of which were computed using five times five-fold cross-validation. The models were compared using the BIC to select the optimal model that fits the data while limiting overfitting. To avoid linear models containing correlated parameters in the same equation, which does not bring much more information than only one and increases the risk of overfitting the data, when two parameters from the subset selected with the hdi function were correlated with r > 0.85, we only kept in the group of tested parameters the one that resulted in the best model. The linear model selection step was applied to select a model for the 60 data points together, as well as for pieces of the cluster mean: the 30 data points of the dry season alone, wet season, irrigated field, rainfed field, as well as the 15 data points of the irrigated field in the dry season, rainfed field in the dry season, irrigated field in the wet season and rainfed field in the wet season.

3. We used these results to form seven piecewise models to be compared to the model selected for the whole dataset (Figure 1B):

   1. model in one piece (calculated with the 60 data points)

   2. dry season + wet season

   3. irrigated field + rainfed field

   4. dry season + irrigated field in the wet season + rainfed field in the wet season

   5. wet season + irrigated field in the dry season + rainfed field in the dry season

   6. irrigated field + rainfed field in the dry season + rainfed field in the wet season

   7. rainfed field + irrigated field in the dry season + irrigated field in the wet season

   8. irrigated field in the dry season + rainfed field in the dry season + irrigated field in the wet season + rainfed field in the wet season

   We chose from these eight models using the BIC. The MSE of composite models was computed by assembling the squared residuals from cross-validation calculated from each distinct linear equation individually into one vector of squared residuals for the whole model and the number of parameters per model was the sum of the number of parameters of each equation.

4. In the case of strong differences in ED responses between fields and seasons, the hdi function might not select ED parameters that can fit each field or season using as output the expression data from both seasons and fields. To take that possibility into account, we repeated the hdi parameter selection with season and field subsets of the data, with B = 300, EV = 2, threshold = 0.7 and fraction = 0.7. We repeated the model selection steps described above using first a pool of the parameters selected using the dry and wet season subsets and then a pool of the parameters selected with the rainfed field subset and the irrigated field subset. We used the BIC to choose from the three models selected from these three groups of parameters.

We used microarray data from two studies (Nagano et al., 2012; Sato et al., 2011) available on the GEO website (accession numbers: GSE36040 and GSE21494), that was already normalized and log-transformed. We used only data from sampling points collected at 8:00, 10:00, 12:00, 14:00 and 16:00 as these times were far enough from sunrise (around 4:30) and sunset (around 19:00). We only used data from samples collected before flowering, from 15 to 62 days after transplanting. Sunrise time varied little enough during that period (from 4:23 to 4:34) that it could be estimated that a series of samples collected at the same time of day would have negligible variation in circadian clock. There were eight samples collected at 8:00, 12 samples at 10:00, eight samples for single replicates and 24 samples for triplicates of eight sampling points at 12:00, eight samples at 14:00 and eight samples at 16:00. We excluded from the analysis genes with an expression value below −7 in more than 17 samples, keeping 19,837 genes, 16,659 of which overlap with the ones in our analysis. We averaged the biological triplicates, resulting in a total of 52 data points. For every gene, we centered individually each of the separate profiles corresponding to a given time of day, thus eliminating potential circadian clock driven differences between these times of sampling. This data constituted the PND.

We tested the transferability of the models determined with the irrigated field data of our experiment, limiting our evaluation to models that were season independent and explained more than half of the variance of the cluster mean. Using the same gene distribution as our irrigated field clusters, we calculated cluster means for the PND. We calculated climatic parameters in the same way as we did for our experiment. The developmental stage parameter was the number of days after transplanting. Input parameters were centered per time of day to match the expression data. We applied the second step of our model selection method to these profiles to select models common to all five time-of-day profiles, using as a subset of parameters the ones in the model determined for our irrigated field data instead of the hdi pre-selection. We only kept a parameter in the new model if it was fitted with a coefficient of the same sign as in the original model.

For the independent analysis of the PND, we used the same method as for our data, which produced 60 gene clusters. The model selection method was a simplified version of the one used for our data. We added to the set of input parameters short-term averages of 8 and 12 hr to account for a possible longer effect of same day conditions, as some samples were collected later in the day than in our experiment. We only selected models common to all five time-of-day profiles, using a unique set of parameters selected from the cluster mean in its entirety and choosing from all possible linear regressions with no more than three parameters using the BIC.

We used the Field Transcriptome Database in Oryza sativa (http://fitdb.dna.affrc.go.jp/) to identify the genes showing a clear environmental response detected by the Nagano et al. models, choosing the ones that had models where both the R² of the overall model and the R²_d of the environmental parameter were over 0.5, with no developmental or circadian terms and an environmental term for either temperature or solar radiation with no gate or a sinusoidal gate and a dose-dependent response.

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

We thank Reynaldo Manuel Jr. for his help with plant growth measurements and collecting the leaf samples, Eloisa Suiton, Julius Borgonia, Godofredo Perez, Josefina Mendoza, Paul Cornelio Maturan, Carlos Casal Jr., and Danilo De Ocampo for their help in the field and the staff of IRRI’s Climate Unit for their help with climate data. We are grateful to Noah Youngs for gathering the GO annotations. We also thank Atsushi Nagano for providing the weather data, detailed result tables and additional details concerning the work published in Naganoet al. (2012). This work was funded by a grant from the National Science Foundation Plant Genome Research Program.

© 2015, Plessis et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.