Effect of phyB and phyC loss-of-function mutations on the wheat transcriptome under short and long day photoperiods

Background Photoperiod signals provide important cues by which plants regulate their growth and development in response to predictable seasonal changes. Phytochromes, a family of red and far-red light receptors, play critical roles in regulating flowering time in response to changing photoperiods. A previous study showed that loss-of-function mutations in either PHYB or PHYC result in large delays in heading time and in the differential regulation of a large number of genes in wheat plants grown in an inductive long day (LD) photoperiod. Results We found that under non-inductive short-day (SD) photoperiods, phyB-null and phyC-null mutants were taller, had a reduced number of tillers, longer and wider leaves, and headed later than wild-type (WT) plants. The delay in heading between WT and phy mutants was greater in LD than in SD, confirming the importance of PHYB and PHYC in accelerating heading date in LDs. Both mutants flowered earlier in SD than LD, the inverse response to that of WT plants. In both SD and LD photoperiods, PHYB regulated more genes than PHYC. We identified subsets of differentially expressed and alternatively spliced genes that were specifically regulated by PHYB and PHYC in either SD or LD photoperiods, and a smaller set of genes that were regulated in both photoperiods. We found that photoperiod had a contrasting effect on transcript levels of the flowering promoting genes VRN-A1 and PPD-B1 in phyB and phyC mutants compared to the WT. Conclusions Our study confirms the major role of both PHYB and PHYC in flowering promotion in LD conditions. Transcriptome characterization revealed an unexpected reversion of the wheat LD plants into SD plants in the phyB-null and phyC-null mutants and identified flowering genes showing significant interactions between phytochromes and photoperiod that may be involved in this phenomenon. Our RNA-seq data provides insight into light signaling pathways in inductive and non-inductive photoperiods and a set of candidate genes to dissect the underlying developmental regulatory networks in wheat.


Background
As sessile organisms, plants must be able to respond to fluctuations in their environment to maximize their reproductive success. To achieve this, plants have evolved a series of regulatory mechanisms to ensure that critical stages of their development coincide with optimal environmental conditions. One important determinant of reproductive success is flowering time, which is strongly influenced by seasonal changes in photoperiod and temperature [1]. In cereal crops, these cues are fundamental to ensure the plant does not flower too early, to prevent exposure of sensitive reproductive tissues to late-spring frosts, or too late, so as to minimize exposure to damaging high temperatures during grain filling [2]. There is a direct link between reproductive success and grain production, so characterizing the regulatory networks underlying flowering time is critical to support the development of resilient crop varieties, to help meet the world's growing demand for food [2].
Plants respond differently to seasonal variation in photoperiod according to the environment to which they are adapted. Whereas some plant species exhibit accelerated flowering in short day photoperiods (SD plants), others flower more rapidly in long days (LD plants). A third class of plants are day-neutral and flower irrespective of the photoperiod. The temperate cereals, including common wheat (Triticum aestivum L.), are LD plants. This ensures that plants remain in a vegetative phase during winter until the lengthening days of spring trigger the irreversible transition to reproductive development [1]. An additional requirement for a long period at low temperatures (vernalization) prevents flowering during the fall, when the days are still relatively long [3].
In wheat and other temperate cereals, the length of the night, rather than the length of the day, is critical for the perception of inductive photoperiods. This has been demonstrated by experiments in which exposing wheat plants to night-breaks (15 m periods of light in the middle of a long night) for at least 12 d was sufficient to accelerate flowering [4]. Loss-of-function mutations in the wheat phytochrome genes PHYTOCHROME B (PHYB) or PHYC, or in the PHOTOPERIOD1 (PPD1) gene abolish the acceleration of flowering by night-breaks, suggesting that these genes are critical to measure the duration of the night [4].
A recent study in Brachypodium proposed a mechanism for the role of these genes in the determination of the photoperiodic response [5]. Phytochromes, a class of red (R,~650 nm) and far-red (FR,~720 nm) light receptors exist as one of two interchangeable forms, P R and P FR . In darkness, the biologically inactive P R form accumulates in the cytoplasm, but upon absorption of R light, P R is converted to the bioactive P FR form and is translocated to the nucleus [6][7][8]. Conversely, exposure to FR light causes the rapid reversion of P FR to the P R form, a reaction that also takes place more gradually during the night (dark or thermal reversion). Therefore, the duration of the night affects the amount of the bioactive P FR form, which has been proposed to be critical for the degradation of the clock protein EARLY FLOW-ERING 3 (ELF3), a direct repressor of PPD1 [5]. High ELF3 protein levels and the repression of PPD1 have been proposed as the main cause of the late flowering phenotypes of the phyC mutant in Brachypodium [5].
PPD1 encodes a PSUEDO-RESPONSE REGULATOR (PRR)-family protein that acts as a positive regulator of flowering in the LD grasses [9][10][11] but as a LD-repressor in the SD grasses rice [12] and sorghum [13], where this gene is referred to as PRR37. In wheat, allelic variation at the PPD1 locus affects photoperiod sensitivity. Whereas the WT Ppd-A1b allele is expressed at very low levels during the night, the Ppd-A1a allele, which carries a promoter deletion encompassing the ELF3 binding site, shows increased transcript levels during the day and, particularly, at night [14]. Wheat varieties that carry the Ppd-A1b allele are referred to as photoperiod sensitive (PS) and those that carry Ppd-A1a as photoperiod insensitive (PI) because they exhibit accelerated heading under SD and reduced differences in heading time between SD and LD. It is important to point out that wheat varieties carrying the PI allele still show a significant acceleration of heading under LD [9,11]. PPD1 induces the expression of FLOW-ERING LOCUS T1 (FT1), which encodes a protein with similarity to the Phosphatidylethanolamine-Binding Protein (PEBP) family [15]. The FT1 protein is translocated through the phloem to the shoot apical meristem, where it forms a hexameric floral activation complex that directly activates the expression of meristem identity genes including VERNALIZATION 1 (VRN1) and FRUITFULL 2 (FUL2). These MADS-box genes play critical roles in triggering reproductive development [16][17][18]. In the cereals, ft1-null mutants exhibit a strong delay in flowering [19].
In addition to their role in the regulation of ELF3, bioactive P FR phytochromes interact in the nucleus directly with PHYTOCHROME INTERACTING FAMILY (PIF) proteins, a class of bHLH transcription factors [20,21]. In Arabidopsis, these interactions induce biochemical changes in the PIF proteins, which result in their ubiquitination and degradation via the 26S proteasome pathway [22]. In this species, accumulating PIF proteins act primarily as negative regulators of light signaling transcriptional networks, so their degradation in response to R light triggers a cascade of photoperiod-mediated transcriptional responses. Despite its important role in the light signaling pathway in Arabidopsis, the role of PIF proteins in the regulation of the photoperiod response in the temperate cereals remains unknown.
Phytochromes can also induce transcriptional variation through modulating alternative splicing (AS) [23]. In Arabidopsis, changes in AS were detected in over 1,500 genes in response to R light, in a PHYB-dependent manner [23]. These target genes include PIF3, whereby greater levels of P FR PHYB increased the frequency of an intron retention event in this gene, disrupting the translated protein's function [24]. In the moss Physcomitrella patens, the phytochrome protein PpPHY4 interacts directly with a splicing regulator to mediate AS in response to light [25]. Previously, the splicing factor RRC was found to mediate phytochrome response in Arabidopsis, suggesting this mechanism may be conserved in angiosperms [26].
Monocot genomes contain three phytochrome genes, PHYA, PHYB and PHYC, with three homeologous copies of each gene in hexaploid wheat [27]. In wheat and Brachypodium, both PHYB and PHYC are required for timely flowering in LD conditions and plants carrying non-functional copies of either phytochrome exhibit extreme delays in flowering, as well as changes in their vegetative morphology [28][29][30]. Using phyB-null and phyC-null Ethyl-methane sulfonate (EMS)-derived mutants in the tetraploid wheat variety 'Kronos', we previously described the sets of genes regulated by PHYB and PHYC in LDs [31]. Despite similar delays in flowering time in both mutants, we found that PHYB regulates approximately six times as many genes as PHYC, and that only a small core of 104 genes were regulated by both phytochromes at the transcriptional level [31]. These commonly regulated genes include several wellcharacterized flowering time genes, such as PPD1 and FT1, and meristem identity genes, including VRN1 and FUL2.
The role of the wheat phytochromes in non-inductive photoperiods remains an open question. Previously, we found that while phyC-null mutants flower later than WT plants in both SD and LD photoperiods, the effect is approximately five-fold smaller in SDs [28]. There is a significant interaction between photoperiod and PHYC, with the WT plants heading earlier in LDs than in SDs, and the phyC-null mutants heading earlier in SDs than LDs [28]. In the current study, we found that phyB-null Kronos mutants also flower significantly earlier in SD than in LD.
To characterize the genes involved in the earlier heading of the phyB-null and phyC-null mutants in SDs than in LDs, we compared the transcriptomes of these mutants under SD and LD conditions. We identified sets of genes regulated by PHYB and PHYC in both SD and LD photoperiods, as well as genes that were regulated only under a specific photoperiod. The findings of this study contribute to our understanding of the complex regulatory networks controlling photoperiod-mediated flowering in wheat.

Results
Effect of phyB-null and phyC-null mutants on heading time We first characterized the effect of Kronos-phyB-null and Kronos-phyC-null mutants on heading time under LD and SD conditions relative to wild-type Kronos (WT), a photoperiod insensitive (Ppd-A1a) spring wheat (Vrn-A1). The WT Kronos headed at 47 d in LD and at 95 d in SD (48 d delay, P < 0.0001), as expected for a LD plant. This result showed that Kronos plants carrying the Ppd-A1a allele still respond to changes in photoperiod.
Under LD, the phyC-null mutant flowered 104 d later than WT and the phyB-null 194 d later than WT ( Fig. 1a-b), confirming the major role of these two genes in promoting flowering under LD [28,31]. By contrast, both phyB-null and phyC-null mutants headed earlier in SD than in LD (108 d earlier for phyB-null, P < 0.001, Fig. 1a, 24 d earlier for phyC-null, P < 0.001, Fig. 1b). This reversal was the result of a much larger delay in heading time in the null mutants under LD (104-195 d) than under SD (31 d and 39 d later than WT, P < 0.0001, Fig. 1a-b). The interactions between photoperiod and genotype were significant for both PHYB and PHYC ( Fig. 1a-b, P < 0.0001).
Consistent with previous characterization of these materials in LDs [28,31], Kronos plants carrying a single null allele in either the A or B homeolog of PHYB or PHYC showed no significant delay in heading date relative to the WT in SDs (Additional file 1, Fig. S1). The same was observed for other traits, so all subsequent results describe comparisons between phyB-null, phyCnull mutants and the WT in a Kronos-PI background.
Effect of phyB-null and phyC-null mutants on plant phenotype under SD We next extended the phenotypic characterization of these mutant lines under SD conditions. Tiller number was significantly lower in both mutants compared to the WT (Fig. 1c), while mean leaf number per tiller was significantly higher in both mutants than in WT plants (Fig. 1d), likely due to the delayed transition of the shoot apical meristem to the reproductive phase. In both phyB-null and phyC-null mutants, flag leaves were significantly longer and wider than WT (Fig. 1e-f).
Stem development was also affected in the phy mutants. Both mutants were significantly taller than WT plants (phyB-null 310 mm taller, P = 9.72 E-06 and phyCnull 220 mm taller, P = 0.00016, Fig. 1g). While the phyB-null and phyC-null mutants did not differ significantly from one another in overall height, their stem structure was markedly different. The phyB-null mutants exhibited a larger number of internodes than either WT (9 more internodes than WT, P = 7.14 E-09 ) or phyC-null mutants (7 more internodes than phyC, P = 3.53 E-07 ), while phyC-null plants had a slightly increased internode number compared to the WT control (2.1 more internodes, P = 0.00013, Fig. 1g). Representative plants of each genotype are shown in Fig. 1h, which was taken when phyB-null mutants reached heading date. Taken together, these results show that both PHYB and PHYC play important roles in regulating vegetative and reproductive development in non-inductive SD conditions.
Characterizing the PHYBand PHYC-regulated wheat transcriptome under SD To investigate the transcriptional changes associated with the later flowering of the phyB-null and phyC-null plants relative to WT in the Kronos background, we performed an RNA-seq experiment in WT, phyB-null and phyC-null plants under SD conditions. We collected tissue from the last fully expanded leaf of four biological replicates per genotype at eight-weeks of age (Additional file 1, Fig. S2). To facilitate comparison with a previous RNA-seq study of the same materials in LD conditions [31], we took samples at the same point of the photoperiod (four hours after dawn). We harvested tissues from eight-week-old plants in our SD experiment so the WT plants were at a similar developmental stage as the WT plants in the LD RNA-seq study, which were sampled at four-weeks of age.
After trimming raw reads for quality and adapter contamination, an average of 45.0 M trimmed 100 bp singleend reads per sample were mapped to unique positions in the IWGSC RefSeq v1.0 genome assembly (Additional file 1, Table S1). Using all normalized read counts mapped to high and low confidence gene models for each sample, we generated a multi-dimensional scaling (MDS) plot (Fig. 2a). Samples grouped into three distinct clusters according to their genotype, reflecting consistent differences in overall transcriptome profile between genotypes and limited differences among biological replicates (Fig. 2a). (c to f) Boxplots represent values of at least five biological replications. Different letters indicate significant differences (Tukey's test P < 0.05). For (a) and (b), * signifies significant differences between photoperiods for each genotype, P < 0.0001. The differences between WT and mutant alleles were also highly significant (P < 0.0001) for both genes and both photoperiods We next performed pairwise comparisons between WT and both mutants to identify PHYB-and PHYCdifferentially expressed (DE) genes under SD conditions. We found that 4.8 times as many genes were regulated by PHYB (7,272 DE genes) than by PHYC (1,511 DE genes, Fig. 2b). Among these DE genes, a greater proportion were positively regulated by PHYB (59.7% with higher expression in WT than phyB-null) than by PHYC (50.6% of genes). There were 815 genes regulated by both PHYB and PHYC, including 783 genes regulated in the same direction and 27 in the opposite direction (upregulated by PHYB and downregulated by PHYC or vice versa, Fig. 2b). Full details of expression data and statistical tests for each pairwise comparison are provided in Additional file 2.
To identify putative functions associated with these transcriptional changes, we performed a GO enrichment analysis for each subset of differentially expressed genes. Among the 7,272 genes regulated by PHYB in SDs, the most significantly enriched terms included 'oxidation-reduction process' and 'protein phosphorylation', while among the 1,511 genes regulated by PHYC, significant terms included 'defense response' and 'cellular iron homeostasis' (Additional file 1, Table S2). In genes commonly regulated by both PHYB and PHYC, enriched terms included 'defense response' and 'protein phosphorylation' (Additional file 1, Table S2).
Changes in development are often associated with differential expression of genes encoding transcription factors. Compared to the overall proportion of genes encoding transcription factors in our dataset (3.2% of 72, 120 expressed genes), an increase was observed for the PHYB-(5.3%) and PHYC-regulated genes (5.4%), and an even larger increase was detected among the genes regulated by both PHYB and PHYC (6.5%). More importantly, several critical genes involved in the regulation of flowering were differentially expressed in phyB and phyC To validate these expression data and to study longerterm trends of the expression of these genes in SD conditions, we performed qRT-PCR analysis for selected candidate genes across six time points, using the same genotypes as in the RNA-seq analysis. At the eight-week time point, the qRT-PCR experiment confirmed the RNA-seq results, showing that transcript levels of VRN1, FT1, FT2, PPD1 and FT3 were all significantly lower in phyB-null and phyC-null mutants compared to WT (Fig. 3). It is important to note that for PPD1, these values represent the combined transcript levels of Ppd-A1a and Ppd-B1b homeologs.
There were also differences in the expression profiles of members of the FT-like family between genotypes. In WT plants, FT1 transcript levels were more than double the levels of FT2 at the 5 w and 8 w time points (Fig. 3), consistent with results from a previous study [32]. By contrast, in both phyB and phyC plants, FT2 was upregulated at an earlier time point (14 w) than FT1, which increased in expression at 17 w in phyC mutants, but remained low throughout the experiment in phyB mutants (Fig. 3). FT3 was expressed at much lower levels than FT1 and FT2, and in the WT both FT-A3 and FT-B3 showed a transient peak in expression at 8 w. In the phyC mutant, FT3 levels started to increase at 14 w and were even higher at 17 w, whereas in the phyB mutant we only observed upregulation of FT-A3 at 17 w (Fig. 3). VRN1 expression increased gradually in both mutant lines throughout this time course, but its transcript levels remained significantly lower than in WT lines at all time points from 5 w onwards (Fig. 3). Taken together, these differences in expression between phytochrome mutants and WT are consistent with the delayed heading date of phyB and phyC mutants compared to WT. These results also confirm that phytochromes play an important role in the regulation of critical flowering genes under both SD and LD photoperiods. The earlier expression of FT2 relative to FT1 and its high transcript levels ( Fig. 3), suggest that this gene may play an important role in promoting wheat flowering time under SD.

Effect of photoperiod on phytochrome-regulated genes
To explore the effect of photoperiod on the differences between WT and the phytochrome mutants, we compared the DE genes generated in the current study in SD collected from 8-week-old plants, with the DE genes in a previous dataset that used the same plant materials grown in LD conditions collected from 4-week-old plants [31]. This SD time point was chosen to match the developmental stages in the WT plants between SD and LD conditions (Waddington stage 3 [33]).
To allow a direct comparison between datasets, we remapped the RNA-seq reads from our earlier LD study to the IWGSC RefSeq v1.0 genome assembly using the same mapping and quantification parameters adjusted for read length. Using this updated genomic reference, 52.8% of all reads mapped uniquely (Additional file 1, Table S3). This LD dataset includes two experimental replicates, each with four biological replications. Genes were considered differentially expressed only when significant in both experiments. This approach reduces the false positive rate, but means that direct comparisons of the number of differentially expressed genes between SD and LD datasets should be approached with caution because SD data represents only a single experimental replicate.
An MDS plot separating all SD and LD samples on the basis of their whole transcriptomic profiles revealed wider differences between photoperiods than between genotypes within a photoperiod (Additional file 1, Fig.  S4). Differences in the expression profiles between experimental replicates in the LD experiment are likely due to variation in the growth chamber light conditions [31]. We identified 3,668 genes that were differentially expressed between WT and phyB-null mutants in both experimental replicates and 424 genes for the corresponding comparisons with the phyC-null mutant. Just 141 of these genes were regulated by both PHYB and PHYC under LD conditions. With slight variations, these results are consistent with our previous study mapping these sequencing data to an older version of the wheat genome [31]. Full details of expression data and An MDS plot with both SD and LD expression data shows the clear differences between the wheat transcriptome of plants grown in different photoperiods (Additional file 1, Fig. S4). In Fig. 4, we divided genes into mutually exclusive classes according to the conditions under which they were differentially expressed between WT and mutant alleles (i.e. regulated by PHYB or PHYC under either SD or LD conditions). For clarity, this figure excludes some pairwise comparisons with low numbers of genes, so the numbers presented in the text do not sum to the complete number of DE genes, which are presented in Additional File 4. In both photoperiods, a greater number of genes were regulated only by PHYB than only by PHYC (Fig. 4). In SDs, 9.6-fold more genes were specifically regulated by PHYB (5,369 genes) than PHYC (561 genes), whereas in LDs, 13.5-fold more genes were specifically regulated by PHYB (2,289 genes) than by PHYC (167 genes, Fig. 4).
There were more genes differentially expressed between WT and mutant genotypes exclusively in SD (589 genes) than exclusively in LD (46 genes, Fig. 4). In addition, the number of genes differentially regulated in a single photoperiod was larger than the number of genes differentially regulated in both photoperiods. For example, there were 1,015 genes regulated by PHYB in both SD and LD, compared to 5,369 and 2,289 genes that were significant in either SD or LD photoperiods, respectively (Fig. 4). Since the LD acceleration of heading time in wheat requires the presence of both PHYB and PHYC, we focused on genes DE in both mutants.
We detected 589 of these DE genes in SD only, 46 in LD only and 43 in both SD and LD (Fig. 4).
In the GO term analysis, significantly enriched functional terms associated with the 43 genes regulated by both phytochromes under SD and LD included 'transcriptional regulation' and 'photoperiodism' (Additional file 1, Table S4). The 24 genes positively regulated by phytochromes (i.e. higher expression in WT than in phy mutants) included FT1, FT2, FT3, PPD-B1, VRN1, FUL2 and FUL3 (Fig. 5, Additional file 1, Table S5). Although the effects were greater in LD, these results confirm that PHYB and PHYC also play a significant role in the activation of these genes in SD in the Kronos-PI background. These results are consistent with our qRT-PCR analysis (Fig. 3). Other genes with the same expression profile as the previous group included a gene encoding a CONSTANS-like CCT-domain protein (TraesCS1A01G220300), and two homeologs encoding MYB-transcription factors with high similarity to RADIALIS (TraesCS6A01G273200 and TraesCS6B01G300600, Fig. 5, Additional file 1, Table S5). One gene (TraesCS1A01G569000LC) was upregulated by PHYB in both SD and LD and by PHYC in SD, but was downregulated by PHYC in LDs (Additional file 4).
Among the 17 genes that were negatively regulated by both PHYB and PHYC in both photoperiods were TraesCS3B01G365300, which encodes a member of the VQ motif protein family of transcriptional regulators, and three genes encoding members of the FLC clade of MADS-box TFs (Fig. 5, Additional file 1, Table S5). Two of these genes encode homeologs of FLC2, which is orthologous to OsMADS51, a SD promoter of flowering in rice [34]. FLC4 encodes the ortholog of ODDSOC2, which functions as a flowering repressor in Brachypodium and is induced by cold treatment in wheat [35]. Interestingly, TraesCS2A01G427200, which encodes WCOR15, a cold responsive gene, was strongly upregulated in both mutant lines, suggesting that phytochromes play an important role in suppressing the cold tolerance pathway in wheat under ambient temperature conditions (Fig. 5, Additional file 1, Table S5). One gene (TraesC-S2A01G019700LC) was downregulated by PHYC in both SD and LD and by PHYC in LD, but was upregulated by PHYB in LD conditions (Additional file 4).
A GO term analysis of the 589 genes regulated by both PHYB and PHYC only under SD, revealed enriched functional terms 'protein phosphorylation' and 'homeostasis' (Additional file 1, Table S4). Among the 367 genes positively regulated within this group, we detected nine WRKY transcription factors, both homeologs of a RADIALIS-like MYB-family transcription factor (TraesC-S7A01G233300 and TraesCS7B01G131600) and TraesCS5B01G054800, which encodes a bHLH TF with similarity to the PIF subfamily (Fig. 6a, Additional file 1, Table S6). We also found in this group FT-A2, FT-A4 and FLC-A1 (Fig. 6a, Additional file 1, Table S6). Among the 213 genes that were negatively regulated by both phytochromes only under SD we identified members of the GATA, G2-like and B-box transcription factor families and TraesCS1A01G334400, which encodes the GA deactivating enzyme GA-2oxidase-A4 (Fig. 6a, Additional file 1, Table S6). The upregulation of four members of the CBF family of cold-activated transcriptional regulators in both phytochrome mutants (Fig. 6a, Additional file 1, Table S6), suggests a similar role to WCOR15 in suppressing the cold tolerance pathway at ambient temperatures, but in this case restricted to SD conditions. Nineteen other genes were either positively regulated by PHYB and negatively regulated by PHYC, or vice versa (Additional file 4). We next studied the 46 genes regulated by both PHYB and PHYC specifically under LD conditions. Among the most significantly enriched functional terms associated with these genes were 'shoot system development', 'long-day photoperiodism' and 'regulation of circadian rhythm' (Additional file 1, Table S4). There were 27 genes positively regulated by both PHYB and PHYC in LD including both homeologs of GIGANTEA, suggesting this gene may play a role in the LD activation of flowering in wheat (Fig. 6b, Additional file 1, Table S7). Among the 16 genes negatively regulated by both phytochromes was TraesCS6B01G315400, which encodes a CONSTANSlike protein, a member of the SPL family of transcription factors and TANDEM ZINC FINGER1, which, in Arabidopsis, interacts with PRR protein components of the circadian clock regulatory network [36] (Fig. 6b, Additional file 1, Table S7). Three other genes were positively regulated by PHYB but negatively regulated by PHYC (Additional file 4).

Effect of genotype on photoperiod regulated genes
Finally, we performed direct pairwise comparisons between SD and LD samples for each genotype (Additional file 5) to identify photoperiod-regulated genes (PRGs). There were a greater number of PRGs in both phyB-null (19,749) and phyC-null (13,740) mutants than in the WT (12,873, Additional file 1, Fig. S5), suggesting that loss-of-function mutations in phyB-null and phyCnull were not sufficient to reduce the large effects on the wheat transcriptome generated by different photoperiods.
Although the different sampling points in LD (4 w) and SD (8 w) were selected so that WT genotypes were at similar developmental stages in both experiments, these results should be interpreted with caution because the effect of photoperiod is conflated with the effect of differences in chronological age. Both phyB-null and phyC-null mutants headed earlier under SD than under LD, so it is likely that the mutant lines were at different stages of development at the time of sampling. This particular sampling strategy likely contributed to the We used this dataset to explore the expression profiles of 19 flowering time genes in different genotypes and photoperiods and the interaction between these factors (Fig. 7).
This analysis confirmed previous results showing that transcript levels of VRN-A1, PPD-B1, FT-B1, FUL-A2, GI, CO-B1 and CO2 are all significantly affected by photoperiod (Fig. 7). Notably, transcript levels of the photoperiod insensitive Ppd-A1a allele were not significantly affected by photoperiod in this dataset, whereas those of the Ppd-B1b allele showed a highly significant effect of photoperiod (P < 0.0001). Transcript levels of GI were more highly expressed in SD, whereas those of CO1 and CO2 were more highly expressed in LD (Fig.  7).
The expression of most of these flowering time genes was also affected by the phyB-null and phyC-null mutations. Significant differences among the three genotypes were accompanied by significant differences between WT and the combined phyB-and phyC-null mutants, with the exception of CO-A2 and CO-B2. The latter result is consistent with a previous study in which CO1 was highly upregulated during the day in phyC-null mutants but CO2 transcript levels were unaffected [28]. The VRN1 paralogs (FUL2 and FUL3) and the florigenrelated genes (FT1 and FT2) all share similar profiles, with higher transcript levels in the WT relative to the phy-null mutants and in LD relative to SD (Fig. 7). FT-A3 transcripts were not detected, whereas FT-B3 transcript levels were higher in SD than in LD, consistent with the known role of this gene as a SD promoter of heading date [19,37].
VRN-A1 and PPD-B1 were the only analyzed flowering promoting genes for which we observed significantly higher transcript levels in the phy-null mutants in SD than in LD and the opposite in the WT (LD > SD). Based on this result, we speculate that these genes could contribute to the earlier flowering of the phy-null mutants in SD than in LD. Expression of these genes was significantly affected by photoperiod and genotype and all three showed significant genotype x photoperiod interactions (Fig. 7). It is important to point out that the SD RNA-seq samples for the phy-null mutants were collected 70-78 days before heading, so they likely represent early stages of flowering induction. It would be interesting to study later time points closer to heading to see if genes that are induced by VRN-A1, such as VRN-B1, FT1, and FT2 [32,38], are upregulated earlier in SD than in LD.

Light signaling and alternative splicing (AS) in wheat
In addition to the differences in transcript levels, we explored whether PHYB or PHYC regulate AS events in wheat using the replicate Multivariate Analysis of Transcript Splicing (rMATS) statistical method [39]. Our RNA-seq datasets show that both PHYB and PHYC regulate the expression of genes encoding components of the splicing machinery (Additional file 1, Table S8). For example, TraesCS2A01G122400, which encodes the large subunit of splicing factor U2AF, was downregulated in phyB-null mutants in both SD and LD conditions and TraesCS1B01G130200, which encodes an Arginine/ serine-rich splicing factor, was upregulated in phyC-null mutants in both SD and LD (Additional file 1, Table S8). There were also several splicing-related genes regulated specifically under SD conditions. Three genes encoding splicing factor subunits were upregulated in both phyB-null and phyC-null mutants, while TraesCS1B01G125800, which encodes pre-mRNA-splicing factor cwc26, was significantly downregulated in both mutants in SD conditions (Additional file 1, Table S8).
To quantify the effect of these changes on AS in wheat, we first identified RNA-seq reads mapping to exon-intron junctions in annotated genes and calculated the frequency of AS events in five different categories (retained intron, skipped exon, alternative 5′ or 3′ splice sites and mutually exclusive exons). Comparing the frequency of each event between WT and mutant genotypes in different photoperiods, we found 5,175 AS events that were significantly regulated by either PHYB or PHYC (FDR P-adj < 0.05). The most commonly observed AS event was intron retention, followed by alternative 3′ splice sites (Fig. 8a).
To classify the events with potentially greater impact on gene function, we looked at the subset annotated genes that showed > 30% variation in their isoform expression levels between genotypes. Among these genes, similar numbers were impacted by AS events in SD and LD (Fig. 8b), although we found a slightly larger number of genes with retained intron events mediated by PHYC in LD conditions (Fig. 8b). Full information on the individual genes impacted by different AS events are provided in Additional file 6.

Discussion
Phytochromes interact with PPD1 in the regulation of wheat heading time Across plant species, one well-characterized function of phytochromes is to regulate flowering time in response to changes in photoperiod. Previous studies have shown that a major function of PHYB and PHYC in the temperate cereals was to accelerate heading time under LD [28][29][30][31], and those results were confirmed here (Fig. 1ab). By contrast, loss-of-function mutations in the orthologous PHY genes in the SD grasses rice and sorghum result in earlier heading under LD [40][41][42][43]. Despite the opposite effect of the phyB-null and phyC-null mutants on heading time in SD-and LD-grasses, these two genes . WT vs. phy indicates an orthogonal contrast comparing the WT versus the two mutants. Data was transformed to provide normality of residuals. **** = P < 0.0001, *** = P < 0.001, ** = P < 0.01, * = P < 0.05, ns = not significant. 1 Since transformation affects the interpretation of the significance of the interactions, we also provide the significance of the interaction in the untransformed data. 2 FT-A3 transcript levels were zero in all samples promote the expression of PPD1/PRR37 in both groups of grasses. The difference between them seems to appear downstream of PHYB and PHYC, since under LD conditions PPD1/PRR37 functions as a flowering repressor in rice [12] and sorghum [13] but as a flowering promoter in the temperate grasses [9][10][11].
One notable result from previous observations on phyC-null mutants [28] and in the current study including both phyB-null and phyC-null mutants is that in both cases the mutants headed earlier in SD than in LD, suggesting that these plants were behaving as if they were SD plants. Since the WT Kronos flowered later in SD than in LD, the differences in heading time between the phy-null and WT were much larger under LD (100-200 d) than under SD (30-40 d) (Fig. 1a-b).
It is important to note that these experiments were all performed in the variety 'Kronos' which carries the PI (Ppd-A1a) allele. This PPD1 allele has a deletion in its promoter region that encompass the binding site of the ELF3 protein repressor [5], resulting in ectopic expression of PPD1 during the night [14]. The night expression of PPD1 is critical for the photoperiodic response, as demonstrated in night-break experiments. Induction of PPD1 in the middle of a 16 h night (SD) by a 15 m pulse of light accelerates heading time almost as much as a LD photoperiod [4]. In Brachypodium, it has been proposed that PHYC activation of PPD1 is mediated by ELF3 [5], so the elimination of an ELF3 binding site in the Ppd-A1a allele in wheat may limit the transmission of the phytochrome signal to PPD1. Therefore, it will be important to determine the effects of phyB-null and phyCnull mutations on heading date in the presence of the PS Ppd-A1b allele to test if the accelerated heading time in SD relative to LD in the phy mutants is maintained in this genetic background. We have initiated the crosses to perform this experiment.
The acceleration of heading time under SD in the phynull mutants has some similarities with SDvernalization, but also some differences. In PS accessions of winter wheat and Brachypodium, an exposure to SD for 6-8 w at room temperature followed by LD replaces the need for vernalization to accelerate heading date [44][45][46][47], but this was not observed in PI wheat accessions [46]. By contrast, we observed SD acceleration in the Kronos-PI background in the presence of phyB-null or phyC-null mutations, which suggests that different regulatory mechanisms are likely involved in these two phenomena.
In order to postulate a mechanism by which heading date is accelerated in SD relative to LD in the phy-null mutants it would be interesting to investigate if the contrasting effect of photoperiod on the regulation of the flowering promoting genes VRN-A1 and PPD-B1 in the phy-null mutants (up regulated in SD) compared to the WT (down regulated in SD) has a role in the earlier flowering of the phy-null mutants in SD relative to LD. The temporal reversion in the order of activation of the FT1 and FT2 genes in the WT and phy mutants may also contribute to the earlier flowering of the phy-null mutants in SD. In the presence of the WT phytochrome alleles, FT1 is expressed to higher levels in Kronos-PI earlier in development than the FT2 gene in SD (Fig. 3) and LD [32]. However, in the phyC-null mutant under SD, FT2 transcripts were upregulated earlier than FT1. By 17 weeks, when these plants were starting to head, FT2 reached very high expression levels (> 10-fold ACTIN) in both the phyB-null and phyC-null mutants. In growth chamber experiments under LD, under SD followed by LD conditions and in fall-planted field experiments, ft2-null alleles conferred only a small delay in heading date [32]. In summary, additional studies will be required to determine the role of PPD1, VRN1 and FT2 in the regulation of heading time under SD in a phy-null background.
Although FT3 transcript levels were lower than other assayed genes, they were also upregulated earlier than FT1 in both phyB-null and phyC-null mutants (Fig. 3). In barley, overexpression of the orthologue HvFT3 accelerates heading in LDs and promotes the transition of the shoot apical meristem from the vegetative to the reproductive stage in both SD and LD [48]. In Brachypodium, BdFTL9, a member of the FT3 clade promotes flowering in SD conditions [47]. This protein forms a floral activation complex only in the absence of BdFT1 (i.e. SD conditions), describing a possible mechanism by which diversity in the PEBP family can finely tune flowering time control according to photoperiod [49]. We identified several other members of the PEBP family that were upregulated in LD conditions (Additional file 5), for which it would be interesting to characterize their role in wheat heading date.
In addition to the PEBP genes, GIGANTEA, VRN2/ GHD7 and CO have been shown to play important roles in the photoperiod response in rice [50,51]. GIGANTEA is a direct promoter of FT in Arabidopsis [52], and in rice GIGANTEA upregulates CO (Hd1) which activates the expression of FT [51,53]. In this study, we show that wheat GIGANTEA was expressed at significantly higher levels under SD than under LD and was positively regulated by both PHYB and PHYC specifically under LD (Fig. 7), suggesting that GIGANTEA may also play a role in the wheat photoperiod pathway. In rice, CO promotes flowering in SD in the presence of functional GHD7/ VRN2 or PRR37/PPD1 alleles, and in LD in the ghd7prr37 double mutant [50] providing an example of how mutations in these photoperiod genes can result in the reversion of the photoperiodic response. Both wheat CO1 homeologs were highly upregulated in both phy mutants, whereas the CO2 homeologs were not affected by the same mutations suggesting that these two paralogs are regulated differently by the phytochrome genes. Interestingly, CO1 transcript levels were higher in SD in the WT and in LD in the phy mutants resulting in a strong interaction between genotype and photoperiod (Fig. 7). We also identified TraesCS7A01G211300, that encodes the ortholog to BdCONSTANS-Like 1 (Additional file 6). This gene is upregulated in LD in WT genotypes, but upregulated in SD in phyB-null and phyCnull mutants. Interestingly this gene was differentially expressed in the Brachypodium elf3-null mutant, suggesting that the TraesCS7A01G211300 and BdCON-STANS-Like 1 orthologs may share similar regulatory mechanisms [5].
We are unable to draw conclusions on the role of the VRN2 locus (duplicated genes ZCCT1 and ZCCT2) because the functional ZCCT-B2a and ZCCT-B2b genes are not annotated in the reference genome used in our study, and the non-functional ZCCT-A1 (TraesCS5A01G541300) and ZCCT-A2 genes (TraesCS5A01G541200) were expressed at low levels (Additional file 5).
When analyzing the expression profiles of flowering time genes it is important to remember that the RNAseq data represent a single time point during the day and during plant development of a very dynamic process of interactions among multiple flowering genes. Therefore, these expression profiles can change if analyzed at different times or developmental stages. Despite this limitation, the information generated for this single time point provided important insights on the complex networks that regulate wheat development in response to the phytochrome signals.

Phytochromes affect plant architecture and vegetative development
In addition to flowering time, we found that mutations in PHYB and PHYC are associated with differences in vegetative development. In both SD and LD photoperiods, the leaves in the phyB-null and phyC-null mutants were longer and wider than in the WT suggesting a more extended or more robust growth (Fig. 1e-f, [28,31]). This is in contrast to the phyC-null mutant phenotype in Brachypodium. The first four leaves of phyC-null plants were shorter than WT in SDs, and not significantly different in length in LD conditions [30]. This discrepancy could be due to the stage of development, since in our study, we measured flag leaves and in Brachypodium, young leaves were studied.
The reduced tiller number in the phyB-null and phyCnull mutants relative to the WT (Fig. 1c) despite their later heading time ( Fig. 1a-b) suggests that PHYB and PHYC may play a role in the promotion of tiller number that is independent of their effect on heading time. A similar reduction in tiller number has been reported in phyB mutants in sorghum [54]. In this mutant, the initiation of axillary meristems and formation of buds occur normally, but bud outgrowth is inhibited compared to the WT from six days after planting. This growth arrest is associated with increased transcript levels of Teosinte Branched1 (SbTB1) and dormancyassociated genes [54,55].
The impact of these alleles on plant height was strikingly different between photoperiods. Whereas under LD conditions both phyB-null and phyC-null mutant lines were shorter than WT plants (51.1 and 59.6 cm shorter, respectively [31]), in SDs both mutants were significantly taller than the WT (31.0 cm for phyB-null and 22.1 cm for phyC-null, Fig. 1g). Interestingly, although overall height of phyB and phyC were similar, the stem development in each mutant was different, with phyBnull mutants exhibiting a greater number of internodes (Fig. 1g). There were several genes regulated by PHYB but not PHYC that may be associated with these phenotypic differences. In both SD and LD, transcript levels of GA20ox-B2 and GA20ox-B4, which encode GA biosynthetic enzymes, were significantly higher in phyB-null mutants than either WT or phyC-null (Additional file 4).
Mutations in the phytochrome genes also affect plant morphology in the short-day grasses. Among rice plants grown in the field under non-inductive LD conditions, those with no functional phytochromes headed earlier, were shorter and had smaller panicles than sister lines with a functional PHYC in a phyA phyB background [42]. In the phyA phyB background, the PhyC gene also affected chlorophyll content, leaf angle and grain size, confirming the multiple pleiotropic effects of the phytochrome mutants in grasses. Sorghum plants carrying non-functional phyB alleles exhibit elongated hypocotyl growth in response to blue light [56] and increased stem elongation and internode number in inductive photoperiods [57]. Similarly, Arabidopsis, a LD-plant, exhibits elongated petioles in phyB mutants [58]. These vegetative phenotypes are characteristic of the shade avoidance response, which is mediated by PHYB in both SD and LD species. The similarities in phenotypes suggest that PHYB may play a role in shade avoidance pathways in wheat that is conserved in the other plant species described above. Some of the multiple genes differentially regulated in PHYB but not in PHYC may play a role in the shade avoidance response.
Our transcriptomic results are also consistent with previous studies that have established a link between phytochromes and the cold regulation pathway [59,60]. In rice, phyB-null mutants exhibit improved cold tolerance [61] and in Arabidopsis, PIF3 binds to the promoters of CBF genes to suppress their expression [62]. We identified four CBF genes and two COR genes that were highly expressed in phytochrome null mutants in SDs (Figs. 5 and 6). Transcript levels of four COR genes were significantly higher in SD than LD in WT and both phy mutants, but the differences were greater in the phy mutants. This demonstrates that in warm ambient temperatures, both PHYB and PHYC act to suppress the activation of the cold responsive pathway during the day. In wheat, a link between light quality and cold tolerance has previously been made [63] and suggests that the destabilization of phytochromes in response to FR light (commonly at higher levels in the dusk) or darkness improves the overall cold tolerance. It would be interesting to test the cold tolerance of the Kronos phytochrome mutants to confirm this activation at the physiological level.

Conclusions
The characterization of loss-of-function mutants for PHYB and PHYC in tetraploid wheat revealed that these genes regulate both vegetative development and heading time, with larger differences between the phy-null mutants and the WT under LD than under SD. Although the major role of PHYB and PHYC in the temperate cereals is the acceleration of heading time under LD, our results also revealed that these genes play a role in repressing flowering under SD. The RNA-seq results presented in this study identified a number of flowering genes showing significant interactions between photoperiod and PHYB/PHYC alleles, which may contribute to the opposite effects of these phytochromes on heading time in different photoperiods. They also revealed subsets of genes that were specifically regulated by PHYB and PHYC in either SD or LD photoperiods, or in both photoperiods, providing insights into light signaling pathways in wheat.

Plant materials, growth conditions and phenotypic measurements
All plant materials for this study were derived from an EMS-mutagenized TILLING population in the tetraploid Triticum turgidum L. ssp. durum Desf. variety 'Kronos' (genomes AABB) which was developed in Dr. Jorge Dubcovsky's lab and described previously [64]. The phyB-null and phyC-null mutants were identified from exome-sequenced lines from this population [65] and were described previously [28,31]. Briefly, we combined null mutations in the A and B homeologs of each gene by marker assisted selection and performed two backcrosses to reduce background mutations. We selfpollinated the mutants for several generations, and used BC 2 F 4 phyB-null and BC 2 F 5 phyC-null mutants for the RNA-seq studies. WT lines correspond to the same Kronos parent used in the backcross. All plants in this experiment carried the Ppd-A1a allele that confers reduced sensitivity to photoperiod [11]. All plants were grown in growth chambers (PGR15, Conviron, Manitoba, Canada) under SD conditions (8 h light/16 h dark) at 20°C day/18°C night temperatures and a light intensity of~260 μM m − 2 s − 1 . All chambers including this SD experiment and a previous LD experiment [31] used similar halide light configurations and were located in the same room.
Heading time was recorded as the number of days after sowing when half of the spike emerged from the boot (Zadoks 55 [66]) using five biological replications (n) per genotype. All physiological measurements were made one time, at maturity. We measured total height and individual internode length (n = 4 for each genotype), total tiller number, leaf number, flag leaf width and length (n = 6 for each genotype). We compared SD data for heading time with previously published heading data of the same mutants grown in the same growth chamber configuration under LD [31].

qRT-PCR assays
Plants used to harvest tissues for qRT-PCR experiments were grown under the conditions described above. Beginning when plants were two-weeks old, we collected tissue from the last fully expanded leaf in liquid nitrogen at three-week intervals until 17 weeks after sowing to cover most of the developmental stages in the mutants. All samples were collected 4 h after the lights came on in the growth chamber (ZT4). We collected four biological replicates of all three genotypes (WT, phyB-null and phyC-null) at each time point, using a single wildtype genotype as a control for both mutants. We extracted RNA using the Spectrum™ Plant Total RNA kit (Sigma-Aldrich, St. Louis, MO) following the manufacturer's instructions. cDNAs were synthesized from 1 μg of total RNA using the High Capacity Reverse Transcription Kit (Applied Biosystems) and quantitative RT-PCR was performed in a 7500 Fast Real-Time PCR system (Applied Biosystems, Foster City, CA) using SYBR Green. Primers for the target genes PPD1 [28], FT1 [15], FT2 [16], FT-A3, FT-B3 [67], VRN1 [15], and the control gene ACTIN [31] were described previously. Expression data are presented as fold-ACTIN levels (molecules of target gene/molecules of ACTIN).

RNA-seq library construction and sequencing
The individual plants used for the RNA-seq experiment were the same plants used for the qRT-PCR and phenotypic studies. For the SD RNA-seq experiment, we extracted RNA samples from four biological replicates of each genotype from eight-week-old plants. One wildtype genotype was used as a control for both phyB-null and phyC-null mutants. At this stage, the apices of the WT plants were at an early stage of spike development (Waddington stage 3 [33]) and the apices of both phyBnull and phyC-null plants were still in the vegetative stage (Waddington stage 1 [33]). Data from the LD RNA-seq experiment was previously described [31] and was generated from RNA extracted from the fullyextended third leaf of four-week-old plants, when the apices of WT plants were at the same developmental stage as in eight-week-old SD-plants. There were two experimental replicates of this data, each comprising 16 samples (four biological replicates of phyB-null, phyCnull and respective wild-type sister lines for each mutant). We assembled and determined the quality of RNA-sequencing libraries using the methods described previously [31]. Libraries were barcoded to allow multiplexing and were sequenced using the 100 bp single read module across two lanes (two biological replicates of each genotype, or six libraries, per lane) on a HiSeq4000 sequencer at the UC Davis Genome Center.

RNA-seq data processing
Raw reads were processed for quality and adapter contamination using the pipeline and parameters described previously [31]. Processed reads were mapped to the IWGSC RefSeq v1.0 genome assembly [68], using GSNAPl [69]. We used parameters -m 4 -n 1 -A sam -N 1 -t 24 for the 100 bp single end read SD data, and parameters -m 2 -n 1 -A sam -N 1 -t 24 for the 50 bp single end read LD data, to generate Sequence Alignment/ Map (SAM) files for each sample. We used high and low confidence gene models from IWGSC Refseq v1.0 annotations. To provide additional context to gene function, we performed a BLASTP search using each annotated gene as a query against the NCBI NR database of proteins. We also added additional annotation information for genes encoding members of different transcription factor families [70], MIKC subclass members of the MADS-box gene family [71] and of the FT-like gene family [67]. Full information of the annotations associated with each differentially expressed gene are provided in Additional file 2.
Uniquely-mapped raw count values were generated using htseq-count (https://github.com/simon-anders/ htseq) as described previously [31]. Genes that showed no raw count values greater than or equal to three in any replicate of any of the three genotypes were discarded, leaving 72,108 genes. The raw counts for these remaining genes were normalized using DESeq2 and differentially expressed genes were detected in pair-wise comparisons between genotypes using a threshold of FDR Padj < 0.01 for both DESeq2 and edgeR statistical tests [72]. For LD data, two experimental replicates were analyzed separately and only genes that were significant in both comparisons (described as "high-confidence" DE genes in our earlier study [31]), were included in this analysis.

Alternative splicing
Alternative splicing events were characterized with rMATS v4.0.1 [39]. A GTF annotation file was created for both SD and LD datasets using Stringtie [73]. Inputs for this file were the sorted BAM files generated during RNA-seq mapping and high and low confidence gene annotations from IWGSC RefSeq v1.1 to specify exonintron boundaries. Genome indices used by rMATS were created from the IWGSC RefSeq v1.0 assembly using STAR (parameter --runMode genomeGenerate) [74]. Fastq files for each sample were trimmed to 100 bp and 50 bp for SD and LD datasets, respectively, using a custom perl script. rMATS was run twice on each dataset, comparing WT with phyB-null and WT with phyC-null samples in both SD and LD datasets, using their respective GTF annotation files [68]. The inclusion level difference for each alternative splicing event was calculated from the number of reads for each replicate that map to a possible inclusion event, normalized by the length of those possible events. The value for each type of event represents the pairwise comparisons of the mean value from four replicates of WT and the respective phy-null genotype. Positive inclusion level differences indicate more reads mapped to an AS event in WT than in the phy-null sample and vice versa. An initial 0.01% splicing difference and FDR < 0.05 filter was used to determine significant alternative splicing events categorized into retained introns, skipped exons, alternative 5′ splice sites, alternative 3′ splice sites, and mutually exclusive exons. A more stringent cutoff of 30% inclusion level difference was used to analyze a subset of these events in greater detail.

Functional annotation
Functional annotation to generate GO terms for each high-confidence and low-confidence gene in the IWGSC RefSeq v1.1 genome was performed as described previously [31].