The transcriptomic response to a short day to long day shift in leaves of the reference legume Medicago truncatula

Growth of plants and tissue sampling

M. truncatula cv ‘Jester’ (Hill, 2000) seeds were scarified by softly scraping them between two pieces of sand paper (grade P160). The seeds were then germinated at 15 °C in gently shaking tubes of water and dark conditions for 24 h. Germinated seeds were then vernalised by being transferred to damp petri dishes and incubated at 4 °C for a further 25 days. The seedlings were subsequently planted in seed raising mix (Daltons Ltd., NZ) in individual cell pots and grown in growth cabinets at 22 °C under 240 µMm⁻²sec⁻¹ cool white fluorescent light. This was in accordance with Institutional Biological Safety Committee approval GMO08-UA006. Soil was kept moist with a complete liquid nutrient media (without Na₂SiO₃; Gibeaut et al., 1997). In the two experiments plants were grown in SD conditions (8 h light/16 h dark) until they were 10 days old and were then shifted at ZT8 into LD conditions (16 h light/8 h dark) for 3 days in experiment 1 (harvest at ZT0 and ZT2) or 5 days in experiment 2 (harvests at ZT4). Three biological replicates were taken each consisting of two pooled trifoliate leaves from different non-adjacent plants giving a total of 18 samples. Only the first trifoliate leaf to unfurl was sampled from a given plant. Samples were immediately frozen in liquid nitrogen.

RNA extraction and sequencing

The 18 frozen leaf samples were ground using five 3 mm sterilised ball bearings per sample and a Geno/Grinder^® 2010 (SPEX^® SamplePrep, USA). Total RNA was then extracted from the ground samples using the RNeasy Plant Mini Kit (Qiagen, Hilden, Germany) following the manufacturer’s instructions and the quantities and qualities of the extracted RNA were then measured using a Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA). Samples were then sent to the Otago Genomics Facility (www.otago.ac.nz/genomics/index.html) and RNA-Seq libraries were prepared. The first experiment (ZT0 and ZT2 samples) used the ScriptSeq Complete Kit (Plant) (100 bp reads; Illumina Inc., USA) while the second experiment (ZT4 samples) used TruSeq Stranded mRNA libraries (120 bp reads; Illumina Inc., USA). Each experiment was sequenced on a single lane of a HiSeq2000 (Illumina Inc., USA).

Complementary DNA synthesis and RT-qPCR analysis

Following RNA extraction, 8 µg of RNA was treated with the TURBO DNA-free Kit (Invitrogen, Foster City, CA, USA) . First-strand complementary DNA (cDNA) was then synthesised using 1 µg of DNase treated RNA using SuperScript IV Reverse Transcriptase (Invitrogen, USA) using an oligo dT primer following the manufacturers instructions. At this point a control reaction where the reverse transcriptase is omitted was run, one reaction per set of replicate samples. When tested alongside synthesised cDNA these reactions control for the presence of genomic DNA contamination.

Measurement of relative abundances of cDNA, as a measure of gene expression, was done using Real time quantitative PCR (RT-qPCR). In this assay 10 µl reactions using Power SYBR Green PCR Master Mix and 2 µl of diluted cDNA (cDNA was diluted 20x prior to RT-qPCR). Control reactions using water, as well as the cDNA reactions which lacked reverse transcriptase, were run simultaneously. RT-qPCR experiments were assembled on 384-well plates and run on an Applied Biosystems 7900HT Sequence Detection System.

Analysis of the RT-qPCR data was done using the 2^−ΔΔCt algorithm (Livak & Schmittgen, 2001). This utilised either Medtr7g089120 a TUBULIN BETA-1 CHAIN gene or Medtr6g084690 a SERINE/THREONINE PROTEIN PHOSPHATASE 2A REGULATORY SUBUNIT (PP2A; previously known as PROTODERMAL FACTOR 2) as reference genes (Kakar et al., 2008). Primers used are available in Table S2.

RNA-Seq analysis

Read trimming of the the FASTQ files was then performed using BBDuk tool in the BBTools suite (v37.54; Bushnell, 2018). This removed the sequencing adapters and low quality sequence (Phred = 20) and retained only reads which were at least 36 bases in length. Furthermore read pairs lacking one of the pair were discarded. Transcripts were then quantified using Salmon (v0.8.2; Patro et al., 2017) which uses quasi-mappings to map the reads to the annotated genes of the Mt4.0v2 transcriptome (Young et al., 2011; Tang et al., 2014). The resulting count tables were then imported into R (R Core Team, 2018) using the tximport package (v1.4.0; Soneson, Love & Robinson, 2015) and principal component analysis (PCA) and DE analysis at the gene level was performed using DESeq2 (v1.16.1; Love, Huber & Anders, 2014). DESeq2 normalizes the data by fitting a negative binomial GLM with a gene-specific dispersion parameter. Clustering was performed using the Mfuzz package (v2.36; Kumar & Futschik, 2007) using minimum centroid distances as a heuristic measure of appropriate cluster number. Briefly this consisted of compromising between cluster size and number by ascertaining when additional clusters only marginally increased the resolution.

Data processing and visualisation

All computation and analysis was done on a Macbook Pro (2012; Intel^® Core^TM i7-3520M upgraded to 16 Gb RAM) from Apple Inc. (Cupertino, CA, USA) running Antergos Linux (v17.12; 4.14.11-1-ARCH kernal). Data processing and visualisation was done in R (R Core Team, 2018) using the Tidyverse suite of packages (v1.2.10; Wickham, 2017) with additional visualisation using the UpsetR package (v1.4; Conway, Lex & Gehlenborg, 2017) and the Superheat package (v0.1; Barter, 2018). Analysis and graphs can be reproduced from the accompanying collection of scripts and files in an accompanying figshare repository (see Thomson, 2018).

Differential expression at each time point

Analysis of this data initially considered pairwise comparisons of gene expression between LD and SD at each timepoint. The significance of any differences observed was assessed using Wald significance tests and since there are three sets of tests, the significance levels of these were adjusted for together using the false discovery rate method.

It was observed that 28,151–29,011 genes had read counts >1 (out of the 50,444 in the Mt4.0v2 transcriptome). Of these 6,824 genes at ZT0 (24% of expressed genes) had statistically different expression (α = 0.05). There were 5,523 genes (19%) at ZT2 and 7,743 genes (25%) at ZT4 (Table S3). When these lists were filtered for those which had >2-fold differences and had >10 mean normalised reads then 2,436, 1,309 and 2,661 genes were judged to be DE at ZT0, ZT2 and ZT4 respectively. These numbers are similar to what was observed in an A. thaliana microarray experiment after a SD to LD shift where 2000 genes were DE (Wigge et al., 2005).

The transcript abundances of the LD induced FT-like genes and five other candidate photoperiod pathway genes were assessed (Fig. 1). While absent in SD, large increases in MtFTa1 transcript abundance were observed in LD at ZT0, ZT2 and ZT4. In addition, appreciable expression of MtFTb1 and MtFTb2 was only seen in LD at ZT4, with minimal to no expression at ZT0 and ZT2 in either SD or LD.

The differences in expression exhibited by the FT-like genes in Figs. 1A to 1C qualitatively agree with previously reported RT-qPCR time course data where in LD MtFTa1 has an approximately constant level of expression across the day (Laurie et al., 2011). MtFTb1 and MtFTb2 have been observed to diurnally peak at ZT4 and ZT16 (Laurie et al., 2011), consistent with the DE observed here at ZT4 (Figs. 1B and 1C). This indicates that despite originating from different experiments, these datasets could be analysed together as a ZT0-to-ZT4 time series to observe the pattern of gene expression change following a SD to LD shift.

To further demonstrate that these datasets could be analysed together we considered the diurnal expression of four other genes measured by RT-qPCR in an independent ZT0-ZT20 time series experiment and compared them to our RNA-Seq data (Fig. 2). Here it was observed that in MtFKF1 (Medtr8g105590) and MtCDF1 (Medtr2g059540) from a low point at ZT0 in SD gene expression increases and peaks at ZT8 (Figs. 2A and 2B). This is also seen in the gene abundances of the RNA-Seq datasets where a large increase in SD at ZT4 is observed. Similar minimal LD expression at these timepoints between the two datasets is also in evidence (Figs. 2E and 2F). Conversely, in the RT-qPCR time course MtCDF2 (Medtr5g041420) and MtCDF4 (Medtr6g012450) have their greatest expression at ZT0 and which in SD sharply decreases at ZT4 (Figs. 2C and 2D) which is also observed in the RNA-Seq transcript abundances (Figs. 2G and 2H).

Overview of the time series analysis

The similarity between the patterns of expression observed in this data with previously reported RT-qPCR results for three genes (Figs. 1A to 1C), as well as independently collected RT-qPCR results for an additional four genes (Fig. 2) indicates that there is no significant batch effect which would bias the interpretation of this data as a time series. We then analysed the data in such a manner.

Across the three timepoints 31,363 genes (out of the 50,444) had a read count >1 at at least one timepoint (see Fig. S1 for dispersion plot and MA-plot) and were included in the analysis. To take an overarching view of the variation between samples, principal component analysis (PCA) was employed with the first two principal components plotted in Fig. 3. Biological replicates clustered together and 64% of the observed variation is explained by the first two principal components which strongly align with the time of sampling (PC1) and the photoperiod condition (PC2).

Typically in a pairwise comparison, differences in transcript abundances can be filtered based on fold-change and levels of expression to focus on genes more likely to be biologically consequential. This is difficult in a time series as it is unknown in this instance which timepoint is most relevant to the regulation of FT-like genes. Reflecting on the gene abundances plotted in Fig. 1 it was observed that the genes with differing expression between LD and SD could be broadly grouped into two classes, those which altered their pattern of expression in response to the photoperiodic shift and those which altered the magnitude of their expression. Specifically, in Fig. 1, the majority of genes (six) change their pattern of expression over time between conditions (MtFTb1, MtFTb2, MtPHYA, MtGI, MtELF3a and MtSOC1a), while two genes maintain a similar pattern, but vary in magnitude of expression between conditions (MtFTa1 and the NF-YC-like gene). It was into these two classes that we decided to segment the time series data followed by clustering.

To classify the genes into these two classes we fit two models to the data. The first model included an interaction term between growth condition (SD or LD) and time of sampling (ZT0, ZT2 and ZT4). This was to identify genes which respond in a condition-specific manner over time such that any changes in the pattern of gene expression brought about by the photoperiodic shift are determined (e.g., MtFTb1 in Fig. 1B). Secondly, we repeated the analysis with a model which lacked this interaction term to test for just the effect of condition (LD vs SD) on the magnitude of gene expression (e.g., the NF-YC-like gene in Fig. 1G). This dual approach facilitated the identification of genes which alter their pattern of expression or only the magnitude of their expression respectively.

Alongside this analysis, an a priori list of 146 candidate genes was also assembled consisting of genes shown to, or are suspected of, participating in flowering time regulation via the photoperiod pathway (Tables 1 and 2). In addition to the FT-like genes, the list incorporates photoreceptors, components of the circadian clock and classes of transcription factor which could potentially regulate the FT-like genes chosen based on their role in A. thaliana. These include candidates from the CDF, TPL and NF-Y families as well as additional transcription factors from families with members known to bind the promoter of FT in A. thaliana compiled by Ridge et al. (2016). These include the genes containing CCT-domain and B-box domains as well as the CIB/BEE-like (CBL), CRYPTOCHROME-INTERACTING BASIC-helix-loop-helix (CIB), PHYTOCHROME INTERACTING FACTOR (PIF), and APETALA2 (AP2) families.

Photoperiod induced changes in the pattern of gene expression

To identify changes in the pattern of gene expression over time, a generalised linear model was fit using the growth condition and time of sampling as predictors along with an interaction term. A likelihood ratio test was then used to test for genes where this interaction term is significant relative to a reduced model lacking the interaction term. A significant result indicates that the expression of the gene responds in a condition-specific manner over time (i.e., a low p-value indicates a change in expression pattern not solely magnitude). This is illustrated in Fig. 1 where MtPHYA, MtELF3a, MtGI, MtFTb1 and MtSOC1a have distinctly different patterns of expression over time in LD compared to SD. On the other hand MtFTa1 and the NF-YC-like gene Medtr1g082660 show a similar pattern in both LD and SD but at differing magnitudes so in this context are not considered to alter their pattern of expression. While MtFTb2 also appears to alter it’s pattern over time, it was not significant in the likelihood ratio test.

This approach resulted in 9,516 genes with altered expression (α = 0.05; full results in Table S4) or 30.34% of those tested with >1 read. To aid interpretation of these changes and quickly identify the timepoint(s) at which individual genes differed in LD relative to SD, Wald significance tests were used to contrast the difference in expression between the two conditions at each of the three timepoints. The significance levels of these genes were adjusted for all three contrasts together using the false discovery rate method. This resulted in 9,427 of the 9,516 genes having differing expression between LD and SD at a minimum of one timepoint with 6,437, 4,511 and 6,159 for ZT0, ZT2 and ZT4 respectively (Fig. S2A). We consider only genes with >2-fold differences with >10 mean normalised reads as DE and there were 3,192 genes meeting this criteria at at least one timepoint. This corresponds to 2131, 1062 and 2168 DE genes at ZT0, ZT2 and ZT4 respectively. Full results are listed in Table S5. In these filtered lists of genes 4.9% of these genes differed in the three timepoints and 32% differed at two or more timepoints demonstrating how most genes have a single peak, predominantly at ZT0 or ZT4 (Fig. 4A). Notably, there were fewer genes DE at ZT2 than the other timepoints and those that did mostly differed at one or both of the other timpeoints too. Only 27% of the 1,062 genes DE at ZT2 are unique to the timepoint (compared to 58.4% for ZT0 and 51.6% for ZT4).

The interaction term was found to be significant in the majority (91/146; 62%) of the candidate genes associated with photoperiodic regulated flowering including 19/24 (79%) of the circadian clock and photoreceptor candidate genes and a striking 8/8 (100%) of the selected AP2 class of genes (Table 1). Of this set of 91 genes, 90 were statistically different (α = 0.05) at one or more timepoints and for 61 genes the difference at at least one timepoint was greater than >2-fold difference with >10 mean normalised reads. For example, a gene encoding a predicted core component of the core circadian clock, MtLHY, is expressed at ∼4-fold higher levels at ZT0 in SD than LD but this situation is reversed at ZT4 when it is ∼7.5-fold higher in LD than SD.

To get a broader view of the results, the mean abundances of all 9,516 genes which altered their pattern of expression were taken and log₂ transformed, before being standardised to have a mean of zero and standard deviation of one. The standardised mean abundances were then clustered into 18 clusters using c-means clustering (Table S6; see Fig. S3A for cluster number optimisation) and visualised in Fig. 4B. This algorithm assigns a membership score (between 0 and 1) for each gene to each cluster describing the degree to which an individual observation belongs to a given cluster (see Fig. S4 for distribution of membership scores). A gene is then assigned to the cluster for which it has the highest membership. Cluster 3, which has 464 genes in it, was of particular interest as MtFTb1 is present. Thus these genes have patterns of expression similar to MtFTb1 and may be involved in regulating MtFTb1 or involved in similar processes. Clusters 9, 11 and 17 (648, 598 and 538 genes respectively) are also of interest as they contain genes with an opposite expression pattern to cluster 3, some of which may thus be negative regulators of MtFTb1 (e.g., MtFTa2 is in cluster 11). Specifically, these three clusters have peaks of expression at ZT4 in SD which are not present in LDs (Fig. 4C).

With regard to candidate photoperiodic flowering time genes, alongside MtFTb1 in cluster 3 there were only two genes from Table 1. These were MtSOC1a (Medtr7g075870) and a BHLH transcription factor gene called MtCBL3 (Medtr7g053410). MtSOC1a is a downstream target of MtFTa1 (Medtr7g084970) and demonstrated to affect flowering time (Jaudal et al., 2018). Another cluster 3 candidate flowering related gene Medtr3g101520 encodes a B3 domain, as does E1, the most important locus in the photoperiod pathway of the tropical legume soybean (Glycine max (L.) Merr.; Xia et al., 2012). On the other hand, the predicted circadian clock-like genes PRR37a, b and PRR59a-c (Matsushika et al., 2000; Nakamichi et al., 2005) all have profiles opposite to that of MtFTb1 with three included in cluster 9 (Table 1). Additionally Medtr1g033620, a SHAQKYF class MYB-like DNA-binding domain gene, is present in cluster 11 which is interesting as proteins of this class have recently been implicated in the regulation of FT in A. thaliana (e.g., EFM and FE; Yan et al., 2014; Abe et al., 2015).

Next, genes in each cluster were then ranked by their fold-change at ZT4 because MtFTb1 has its greatest difference in expression between the two conditions at this timepoint (962-fold up in LD; Fig. 1B). We focused on transcription factor genes and observed significant fold changes in transcript levels of genes encoding a number of zinc finger proteins. For instance, in cluster 3 at ZT4, the zinc finger (Ran-binding) family gene Medtr4g113840 was 45-fold elevated in LD compared to SD and the DOF-type zinc finger DNA-binding family gene Medtr3g091820 was increased 2-fold. In contrast, in cluster 9 the B-box type zinc finger protein Medtr2g073370 was 29-fold more abundant in SD than in LD at ZT4 and in cluster 11, Medtr2g059540, also a DOF domain zinc finger gene (denoted MtCDF1 in Table 1) was ∼300-fold higher in SD than in LD. Furthermore in cluster 17, the zinc finger (Ran-binding) family gene Medtr6g069400 is 7-fold elevated in SD compared to LD at ZT4.

Changes in the magnitude of gene expression

There is a class of genes which change their level of expression in response to the shift from SD to LD conditions, but this does not alter the relative changes which occur at differing timepoints. Thus it is only the magnitude, not the pattern of expression which changes (e.g., MtFTa1 or the NF-YC-like gene Medtr1g082660 in Figs. 1A and 1G). To identify these genes a simpler model was fit the data which lacked the interaction term between growth condition and time of sampling.

It was observed that in this model 8,695 genes differed in the magnitude of their expression between conditions (α = 0.05; Table S7), but this was reduced to 4,694 when those that also altered the pattern of their expression were omitted. Therefore 14.96% of genes with >1 read (4,694/31,363 genes) alter just the magnitude of their expression in response to the shift in photoperiod conditions. To provide timepoint level resolution of these changes Wald significance tests were again used to contrast the 4,694 genes (Table S8) with the significance levels of these genes adjusted for all three contrasts together using the false discovery rate method. This resulted in 4,161 of the 4,694 genes differing in magnitude at least one timepoint with 2,715, 2,268 and 2,666 for ZT0, ZT2 and ZT4 respectively (Fig. S2B). When considering genes with >2 fold differences and >10 mean normalised reads as DE, only 819 genes differed in magnitude at at least one timepoint and at ZT0, ZT2 and ZT4 this corresponded to 457, 353 and 454 genes respectively. Here it was observed that only 65 genes (7.9%) of this set of genes are consistently higher in LDs than SDs while 54 genes (9.99%) are consistently lower (Fig. 5A). The numbers in these classes go up to 187 genes (22.8%) and 139 genes (17%) respectively when genes differing at two or more timepoints are included. Like the class of genes which altered their pattern of expression there are fewer DE genes at ZT2.

Results for the candidate photoperiod loci are summarised in Table 2. Specifically these are the candidate photoperiodic flowering time genes which were not classed as changing their pattern of expression between conditions. Here 15/55 (27%) of the genes are classed as altering just their magnitude in response to the photoperiodic shift. However, only 9/15 (60%) have statistically different levels of expression at two or more timepoints and looking down the list MtFTa1 is the only gene to consistently differ >2-fold with >10 mean normalised reads. Thus, none of these candidate photoperiod genes are expressed similarly to MtFTa1 which shows >30-fold higher levels in LD than in SD at all timepoints.

We then assessed the sets of genes whose expression is either consistently higher or lower in LDs than SDs (65 and 54 genes respectively) for other candidate genes which could potentially play a role in the regulation of flowering time. These include Medtr8g091720 a NF-Y-like gene, which exhibits consistently greater expression (2.5–5.6-fold) in LD than in SD, like that of MtFTa1. Genes which consistently show reduced expression in LD compared to SD include Medtr2g014200 which encodes a squamosa promoter-binding-like protein as well as genes associated with sugar transport. For instance, Medtr3g074180 encoding a trehalose-6-phosphate phosphatase is 2–3 fold lower accross the three timepoints and Medtr0204s0040 a sugar porter family MFS transporter which not expressed in LD at all. Sugar transport is a process linked to flowering time and the regulation of FT in A. thaliana (Wahl et al., 2013).

In addition, while not previously linked to flowering time, these lists also contain a number of genes which likely have regulatory functions such as Medtr5g079220 encoding a R2R3-MYB transcription factor, Medtr3g107940 which produces a FBD protein and Medtr8g012655 which is the gene for an ethylene response factor. These genes all have greater expression in LD compared to SD. Conversely, genes which have reduced expression in LD compared to SD include Medtr7g012790 encoding a circadian clock coupling factor ZGT, Medtr7g105780 encoding an ovate transcriptional repressor, Medtr3g031220 which produces a WRKY transcription factor and Medtr8g026960 which is the gene for a homeobox associated leucine zipper protein.

Relaxing the criteria for differential expression slightly and reconsidering the candidate photoperiod loci also suggests that MtCDFf (Medtr6g027450) could be investigated further as it is higher in LD in all three timepoints (8.57-fold, 3.48-fold and 12.99-fold respectively) however its expression is overall quite low with an average of only 2.7 normalised reads at ZT4. Other potential candidates on the list include the ELF4-like gene Medtr8g020200, consistently 1.8–2.3-fold higher in SD compared to LD, as is MtCMF17 (Medtr1g044785) although, like MtCDFf, the expression in these datasets is low. This may reflect that these genes are cell-type specific and so only expressed in a fraction of those sampled.

All 4,694 genes which differed in the magnitude of their expression between conditions were then clustered into 18 clusters (Fig. 5B, Fig. S3B and Table S9). MtFTa1 was present in cluster 1 which had 331 genes. However, none of the selected photoperiodic candidate genes clustered with MtFTa1. Clusters 9 (259 genes) and 10 (239 genes) have patterns opposite to that of cluster 1 (Fig. 5C). A TPL-like gene from Table 2 is in cluster 9 and a second TPL-like gene is in cluster 10. In addition, in clusters 9 and 10, a number of other flowering time candidates are present including Medtr0020s0120 in cluster 9, which is similar to the FT antagonist TERMINAL FLOWER 1 in A. thaliana (Jaeger et al., 2013; Wickland & Hanzawa, 2015). Also present in cluster 9 are a trio of genes encoding B3 binding domain proteins Medtr1g021410, Medtr1g021435 and Medtr1g021500 and pair of genes which encode SHAQKYF class MYB transcription factors Medtr0036s0260 and Medtr5g027550. These genes are notable for genes containing these domains have been associated with flowering time, as has the jumonji domain protein encoding gene Medtr2g011630 in cluster 10 (Xia et al., 2012; Abe et al., 2015; Yan et al., 2014).

A total of 35/65 of the consistently differentially expressed genes which are higher in LDs than SDs (Fig. 5A) are present in cluster 1. These include Medtr1g099440 which encodes a membrane-associated kinase regulator-like protein, Medtr6g086805 a heat shock transcription factor gene and Medtr4g009110 which encodes a helix loop helix DNA-binding domain protein. Similarly, clusters 9 and 10 contain 19/54 of the consistently differentially expressed genes which are lower in LDs than SDs. These genes include the ethylene response factor gene Medtr5g016750 and Medtr4g119422 encoding a cullin-like protein.

We then ranked the clustered genes by their fold-change between LD and SD at ZT0. Strikingly, in cluster 1, within the top 13 genes with the largest fold increase in LD compared to SD at ZT0, MtFTa1 is ranked 8th, while 10 of the other genes were homologues of IRON MAN (IMA)/FE-UPTAKE-INDUCING PEPTIDE 1 (FEP1) genes. These encode mobile signalling peptides integral for the uptake of of iron from the soil and octuple ima/fep1 mutants in A. thaliana result in severe chlorosis (Grillet et al., 2018; Hirayama et al., 2018). They were also identified in a recent reannotation of the M. truncatula genome to identify small, secreted peptides (De Bang et al., 2017). In addition cluster 3, which is similar to cluster 1, has the remaining annotated IMA/FEP1 genes as 4 of the top 5 genes with the largest fold-change differences at ZT0.