Conservation and Diversification of Circadian Rhythmicity Between a Model Crassulacean Acid Metabolism Plant Kalanchoë fedtschenkoi and a Model C3 Photosynthesis Plant Arabidopsis thaliana

Crassulacean acid metabolism (CAM) improves photosynthetic efficiency under limited water availability relative to C3 photosynthesis. It is widely accepted that CAM plants have evolved from C3 plants and it is hypothesized that CAM is under the control of the internal circadian clock. However, the role that the circadian clock plays in the evolution of CAM is not well understood. To identify the molecular basis of circadian control over CAM evolution, rhythmic gene sets were identified in a CAM model plant species (Kalanchoë fedtschenkoi) and a C3 model plant species (Arabidopsis thaliana) through analysis of diel time-course gene expression data using multiple periodicity detection algorithms. Based on protein sequences, ortholog groups were constructed containing genes from each of these two species. The ortholog groups were categorized into five gene sets based on conservation and diversification of rhythmic gene expression. Interestingly, minimal functional overlap was observed when comparing the rhythmic gene sets of each species. Specifcally, metabolic processes were enriched in the gene set under circadian control in K. fedtschenkoi and numerous genes were found to have retained or gained rhythmic expression in K. fedtsechenkoi. Additonally, several rhythmic orthologs, including CAM-related orthologs, displayed phase shifts between species. Results of this analysis point to several mechanisms by which the circadian clock plays a role in the evolution of CAM. These genes provide a set of testable hypotheses for future experiments.

Crassulacean acid metabolism (CAM) improves photosynthetic efficiency under limited water availability relative to C 3 photosynthesis. It is widely accepted that CAM plants have evolved from C 3 plants and it is hypothesized that CAM is under the control of the internal circadian clock. However, the role that the circadian clock plays in the evolution of CAM is not well understood. To identify the molecular basis of circadian control over CAM evolution, rhythmic gene sets were identified in a CAM model plant species (Kalanchoë fedtschenkoi) and a C 3 model plant species (Arabidopsis thaliana) through analysis of diel time-course gene expression data using multiple periodicity detection algorithms. Based on protein sequences, ortholog groups were constructed containing genes from each of these two species. The ortholog groups were categorized into five gene sets based on conservation and diversification of rhythmic gene expression. Interestingly, minimal functional overlap was observed when comparing the rhythmic gene sets of each species. Specifcally, metabolic processes were enriched in the gene set under circadian control in K. fedtschenkoi and numerous genes were found to have retained or gained rhythmic expression in K. fedtsechenkoi. Additonally, several rhythmic orthologs, including CAM-related orthologs, displayed phase shifts between species. Results of this analysis point to several mechanisms by which the circadian clock plays a role in the evolution of CAM. These genes provide a set of testable hypotheses for future experiments.

INTRODUCTION
All organisms confront environmental fluctuations on a daily basis, including changes in light, temperature, predation risk and water availability. Many of these changes exhibit predictable diurnal oscillations, which incentivized organisms to evolve internal timekeeping systems in anticipation of environmental changes. As a result, an emergent circadian clock enabled organisms to synchronize their biological activity with external cues (Bell-Pederson et al., 2005). The circadian clock is important for plant survival, due to their sessile habit, and controls time-of-day specific biology. Greenham and McClung (2015) demonstrated that defense pathways become subordinate to the circadian clock, allowing the plant to activate defense pathways to the time of day when the threat posed by herbivores is maximal. Numerous processes have been documented to be under the control of the circadian clock; however, the extent to which the clock controls basic physiology is still not well understood (Piechulla, 1988;Merida et al., 1999;Thain et al., 2000;Hotta et al., 2007;Mallona et al., 2011;McClung, 2013;Ming et al., 2015).
Crassulacean acid metabolism (CAM) provides a unique physiological system for the study of circadian control over physiology due to a temporal separation of CO 2 fixation and inverted day/night stomatal movement patterns. It is widely accepted that CAM photosynthesis plants have evolved from C 3 photosynthesis plants (Yang et al., 2015Yin et al., 2018). CAM photosynthesis results in improved photosynthetic and water-use efficiency, relative to C 3 photosynthesis, and it is hypothesized that the strict temporal control of CAM processes is maintained by the circadian clock in CAM plants (Wyka et al., 2004;Boxall et al., 2005;Hartwell, 2005Hartwell, , 2006von Caemmerer and Griffiths, 2009;Ming et al., 2015). Evidence from several levels of biology have been reported supporting this hypothesis. At the physiological level, Wilkins (1992) demonstrated that stomatal conductance in CAM plants displays robust 24 h rhythms in constant light and temperature, and at the metabolic level, CO 2 uptake and internal CO 2 concentrations displayed robust rhythms in similar conditions. Hartwell et al. (1999) reported circadian rhythms at the molecular level by showing 24 h rhythms in transcript abundances of one key CAM gene phosphoenolpyruvate carboxylase kinase. Furthermore, at the biochemical level, phosphorylation of another key CAM gene, phosphoenolpyruvate carboxylase, by phosphoenolpyruvate carboxylase kinase displays circadian rhythms in constant darkness (Nimmo et al., 1984). However, a thorough genomic investigation of circadian control in CAM plants has not been performed.
Identification and characterization of circadian-clock regulated genes requires measuring gene expression over a 24-h period (Bar-Joseph et al., 2012;Li et al., 2015;Hughes et al., 2017). Moreover, a large amount of time-course data, with 2to 4-h sampling intervals, are currently available: Arabidopsis thaliana (Mockler et al., 2007), Oryza sativa (Filichkin et al., 2011), Populus trichocarpa (Filichkin et al., 2011), Solanum lycopersicum (Higashi et al., 2016), Brachypodium distachyon (Koda et al., 2017), and Ananas comosus (Ming et al., 2015). Furthermore, numerous algorithms exist for detecting rhythmic expression within such data sets; however, different algorithms yield inconsistent results when used on a common data set (Keegan et al., 2007;Doherty and Kay, 2010;Deckard et al., 2013). For instance, the definition of "periodic" can vary between different algorithms, and depending on the algorithm used, can result in different genes being identified as rhythmic. Additionally, distributions of scores from periodicity detection algorithms are not bi-model, such that there is no clear distinction between periodic and non-periodic genes. Recently, a new approach for detection of rhythmic gene expression in time-course data was developed for integrating multiple rhythmic detection algorithms (Kelliher et al., 2016), allowing one to leverage multiple definitions of "periodic" as detailed by the different algorithms and arrive at "high confidence" rhythmic gene sets that don't necessarily represent the only rhythmic genes. Applying such a method to C 3 and CAM photosynthesis plants will provide insights into how the circadian clock propagates environmental signals across the two different photosynthetic types as well as provide clues to how CAM evolved from C 3 .
Therefore, the aim of this study is to investigate the conservation and diversification of rhythmic gene expression between a model C 3 plant A. thaliana and a model obligate CAM plant K. fedtschenkoi. Through comparative analysis of time-course data using rhythmic detection algorithms (Deckard et al., 2013) we identified commonalities and differences in rhythmic genes between A. thaliana and K. fedtschenkoi. Our results provide further lines of evidence supporting the hypothesis that the circadian clock played a role in the evolution of CAM.

Time-Course Gene Expression Data
The diel expression data for Kalanchoë fedtschenkoi and Arabidopsis thaliana were obtained from Yang et al. (2017) and Mockler et al. (2007), respectively. K. fedtschenkoi genes with a max FPKM < 1 were removed. The K. fedtschenkoi expression data were collected at 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, and 24 h whereas the A. thaliana data were collected at 4, 8, 12, 16, 20, and 24 h after the initiation of the light period. Since the A. thaliana gene expression data were measured at 4-h intervals and the K. fedtschenkoi data were collected at 2-h intervals, the A. thaliana data were imputed to arrive at expression profiles for all A. thaliana and K. fedtschenkoi genes on the same time scale. Interpolation of the A. thaliana data for comparative gene expression analysis has been conducted before . In the present study, the piecewise cubic Hermite interpolating polynomial (pchip) interpolation function in the pandas Python library was used to interpolate the A. thaliana data at time points consistent with time intervals : 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, and 24 h. Pchip was used to maintain the shape of the data (Dong et al., 2011;Luo et al., 2014).

Identify Rhythmic Genes in K. fedtschenkoi and A. thaliana
The rhythmic genes in K. fedtschenkoi and A. thaliana were identified by using a methodology developed by Kelliher et al. (2016), which applied four periodicity detection algorithms, i.e., de Lichtenberg (de Lichtenberg et al., 2005), JTK-CYCLE (JTK) (Hughes et al., 2010), Lomb-Scargle (Scargle, 1982) and persistent homology (Lomb, 1975), to the time-course gene expression data. The initial step in the overall method is to apply the periodicity detection algorithms to each dataset, summing the gene rankings output from each algorithm and reranking based on the cumulative ranks. We did not use the persistent homology step as it has been shown to give poor results when applied to higher organisms (Deckard et al., 2013). The second step in the Kelliher et al. (2016) method is to visually assess rhythmicity in the dataset of cumulative ranked genes. A 500 gene cohort was plotted to visually inspect where rhythmicity in expression decreased. A step size of 500 genes was empirically determined. A cutoff of 10,000 and 10,500 was used to identify top ranked genes for A. thaliana and K. fedtschenkoi, respectively. To validate this approach on time-course gene expression data with the length of one period, the top 10,000 ranked genes in A. thaliana were compared to a previously published rhythmic gene set using the same data and found a 48% overlap (Mockler et al., 2007). A 100% overlap was not expected because several factors such as data normalization, periodicity detection and/or algorithm score cutoff can result in varying sets of genes even from the same data set (Keegan et al., 2007;Kelliher et al., 2016).
We then applied a threshold based on the JTK algorithm at a p ≤ 0.06. Kelliher et al. (2016) based a threshold on the Lomb-Scargle algorithm but the JTK algorithm has been shown to give preference to cosine and peaked gene expression profiles (Deckard et al., 2013) and a previous study examining rhythmicity in transcription factors and transcription coregulators in the CAM plant pineapple found cosine and peaked profiles were most successful in identifying cycling transcription factors and transcription co-regulators (Sharma et al., 2017).

Functional Analysis of Rhythmic Genes K. fedtschenkoi and A. thaliana
Gene Ontology (GO) terms for the K. fedtschenkoi and A. thaliana genes were used from the gene annotation information downloaded from Phytozome v12.1 (Goodstein et al., 2012). K. fedtschenkoi genes encoding putative transcription factors were retrieved from Yang et al. (2017), while A. thaliana transcription factors (TFs) were retrieved from PlanTFDB v4.0 (Jin et al., 2017).
Using ClueGO (Bindea et al., 2009), observed GO biological process were subjected to the right-sided hypergeometric enrichment test at medium network specificity selection and p-value correction was performed using the Holm-Bonferroni step-down method (Holm, 1979). There was a minimum of 3 and a maximum of 8 selected GO-tree levels and each cluster was prescribed to include a minimum of between 3 and 4% of genes associated with each term. GO-term fusion and grouping settings were selected to minimize GO-term redundancy and the term enriched at the highest level of significance was used as the representative term for each functional cluster. The GO terms with p-values less than or equal to 0.05 were considered significantly enriched.

Identification of Orthologous Genes
Ortholog groups consisting of proteins from 26 plant species were obtained from Yang et al. (2017). In short, ortholog groups were constructed using FastOrtho with default parameters selected, except for a BLASTp E-value cutoff of 1e-5 and an inflation value of 1.3. To determine if K. fedtschenkoi genes have homologs in other plant species, K. fedtschenkoi protein sequences were queried a local plant proteome database using BLASTp, with an E-value cutoff of 1e-10, implemented in the National Center for Biotechnology Information (NCBI) BLAST+ (Camacho et al., 2009). This local proteome database contained protein sequences from 83 plant species (Supplementary Table S3), including 70 species from PLAZA 4.0 (Van Bel et al., 2017) and 10 species from Phytozome 12.0 (Goodstein et al., 2012) and three Agave species with transcriptome sequenced (Gross et al., 2013;Abraham et al., 2016).

Identifying Rhythmic Genes in A. thaliana and K. fedtschenkoi
Using the Kelliher et al. (2016) method to identify "high confidence" rhythmic gene sets, 9,338 genes in A. thaliana ( Figure 1A) and 8,769 genes in K. fedtschenkoi ( Figure 1B) were estimated to be rhythmically expressed with a 24-h period, representing 34% of the genes in the A. thaliana genome and 28% of the genes in the K. fedtschenkoi genome. The expression patterns in Figures 1A,B reveal a continuum of phase-specific gene expression, which is a typical behavior of periodic genes controlled by the circadian clock in A. thaliana (Covington et al., 2008), Solanum lycopersicum (Higashi et al., 2016) and Brachypodium distachyon (Koda et al., 2017). Multiple criteria for evaluating expression patterns for individual genes in A. thaliana and K. fedtschenkoi in our study are provided in Supplementary Tables S1, S2.

Functional Comparison of Rhythmic Genes Between A. thaliana and K. fedtschenkoi
A functional enrichment of both rhythmic gene sets was performed using ClueGO (Bindea et al., 2009) to gain insight into the physiological relevance of the circadian clock. Top functional groups enriched in the A. thaliana rhythmic gene set include pigment biosynthetic processes, ribonucleoprotein complex biogenesis and response to high light intensity, red or far red light, light stimulus and cold (Supplementary Table S4). The K. fedtschenkoi rhythmic gene set was enriched for organic acid metabolic processes, carboxylic acid metabolic processes, protein targeting to chloroplast, starch metabolic processes, oxylipin metabolic processes, phosphate-containing compound metabolic FIGURE 1 | Rhythmically expressed genes in Arabidopsis thaliana and Kalanchoë fedtschenkoi. Taking the intersection of a visual inspection of cumulatively ranked genes and a JTK-CYCLE cutoff of 0.06 resulted in a qualitative selection of highly rhythmic genes in each species. (A) 9,338 genes were identified as highly rhythmic in A. thaliana indicated by the wave of peak expression through time. (B) 8,769 genes were identified as highly rhythmic in Kalanchoë fedtschenkoi as indicated by the wave of peak expression though time. White and black bars indicate daytime (12-h) and nighttime (12-h), respectively. The y-axes represent individual genes and transcript levels are depicted as a z-score change relative to mean expression for each gene, where values represent the number of standard deviations away from the mean. (C) Types of rhythmic orthologs and genes in A. thaliana (left) and K. fedtschenkoi (right). Type 1 refers to orthologs that are rhythmic in both species. Type 2 refers to orthologs that are rhythmic only in A. thaliana. Type 3 refers to orthologs that are rhythmic only in K. fedtschenkoi. Type 4 refers to genes found to be rhythmic only in A. thaliana. Type 5 refers to genes found to be rhythmic only in K. fedtschenkoi. Types are designated by their type names and separated by red dotted lines. Vertical numbers are the number of genes in each group from each plant species. processes, photosynthesis (light harvesting), glyceraldehyde-3phosphate metabolic process, isoprenoid biosynthetic process and plastid organization ( Supplementary Table S4). Strikingly, no overlap in enriched functional groups was observed when comparing the dominate rhythmic processes between these two species.
To further investigate the physiological relevance of the circadian clock in a temporal context, a chronological gene co-expression network (GCN) was constructed from each rhythmic dataset, with subnetworks representing each phase of the circadian system (Supplementary Figure S2). Spearman's rank correlation coefficient was calculated for all pair-wise combinations of the rhythmic genes in each dataset. The A. thaliana and K. fedtschenkoi GCN contained 9,338 rhythmic genes linked by 4,174,924 edges and 8,769 rhythmic genes linked by 4,590,438 edges, respectively, with a r ≥ 0.8 as a cutoff of co-expression. Visualizing phase calls of every gene with a color scheme, rhythmic genes can be seen chronologically connected in a circular network (Figure 2A), consistent with the diagonal stripe observed in the rhythmic gene expression heat maps for each species (Figures 1A,B). A similar result has been reported in a previous study in B. distachyon (Koda et al., 2017).
In the rhythmic gene sets for A. thaliana and K. fedtschenkoi, 644 and 923 genes were found to encode putative transcription factors (TFs), respectively. Consistent with the previous study (Koda et al., 2017), we found that the rhythmic expression of genes is correlated to transcription factors that were also rhythmically expressed at the transcript level. Rhythmic TFs in both species were found in all phases of the day and night (Supplementary Figure S2) and only photosynthesis and starch metabolism were common to both species' rhythmic gene lists (Supplementary Table S5).
To evaluate the conservation of rhythmicity and phase of expression between A. thaliana and K. fedtschenkoi, ortholog groups were constructed using the protein sequences from these two species (Figure 1C). Between the rhythmic gene sets, five types of orthologs and genes were identified: Type 1 containing orthologs that are rhythmic in both species, Type 2 containing orthologs that are rhythmic only in A. thaliana, Type 3 containing orthologs that are rhythmic only in K. fedtschenkoi, Type 4 containing genes that are rhythmic and found only in A. thaliana, and Type 5 containing genes that are rhythmic and found only in K. fedtschenkoi. Gene lists for each species in each ortholog group type can be found in Supplementary Table S6. Gene numbers and functional investigations in each type are described below.
Type 1: Rhythmic Ortholog Groups Shared Between A. thaliana and K. fedtschenkoi 2,641 ortholog groups were identified to contain at least one rhythmic gene in each of the two plant species analyzed here. GO enrichment of each subset of Type 1 rhythmic genes displayed a 37% overlap in enriched biological processes, including protein phosphorylation, organic acid metabolic process, transport, photosynthesis and generation of precursor metabolites and energy (Supplementary Table S7). In general, cross-species identification of rhythmic orthologs can be challenging due to "many-to-many" orthologous relationships and when an ortholog group contains genes which are rhythmically expressed in only one of the two species. Therefore, further analysis of conserved rhythmicity within the Type 1 group was limited to only those ortholog groups displaying a single expression profile within a species. This restriction resulted in 768 ortholog groups between the two species where both species' genes were rhythmic. Rhythmic genes were observed in all phases of the day and night that were sampled in each species (Figure 2A). Several key CAM genes involved in CO 2 fixation in K. fedtschenkoi  were found to be rhythmic, as well as several of their A. thaliana orthologs (Supplementary Tables S8, S9). Evidence of both synchronization and desynchronization were found in one-to-one ortholog groups, e.g., with the NADP-malic enzyme (Supplementary Table S9).
To further understand the functional relevance of the rhythmic genes, phase shifts between rhythmic orthologs in the two species were examined next. Phase shifts are defined here as the difference between two rhythmic genes' phase calls, which is the time point where maximum transcript abundance occurs for a gene. We define two orthologs to be synchronized if their phase calls are within 4 h of each other and unsynchronized otherwise. Of the 768 ortholog groups, 386 (50%) ortholog groups contained ortholog pairs with synchronized gene expression ( Figure 2B), with several synchronized ortholog pairs phased between morning and evening ( Figure 2C). The second sub-type, unsynchronized orthologous gene pairs, contained 382 (50%) ortholog groups ( Figure 2B). One group of K. fedtschenkoi genes peaking at midday had orthologs in A. thaliana peaking in the morning, while another group of K. fedtschenkoi genes peaking at midday had orthologs in A. thaliana peaking just after dusk, suggesting differential regulation of rhythmic orthologous gene pairs between A. thaliana and K. fedtschenkoi. TFs were also found in the 768 ortholog groups and in all phases of the day and night in each species (Figure 2D). Of the 61 TF ortholog groups, 40 (66%) contained unsynchronized TF ortholog pairs ( Figure 2E). The 61 orthologous TFs in A. thaliana and K. fedtschenkoi are listed in Table 1, along with the respective A. thaliana gene annotation. Several orthologous TFs had similar phase shifts as the two groups of differential regulated orthologs described earlier (Figure 2F), potentially accounting for the difference in gene expression between orthologous gene pairs. Type 2: Rhythmic A. thaliana Genes With Arrhythmic Orthologs in K. fedtschenkoi 2,797 ortholog groups were identified that had at least one A. thaliana gene that was rhythmic and all K. fedtschenkoi genes were arrhythmic. A total of 3,275 A. thaliana genes found to be rhythmic that had no rhythmic K. fedtschenkoi orthologs ( Figure 1C). Only two biological processes were enriched in this A. thaliana gene set: cellular nitrogen compound metabolism and response to endoplasmic reticulum stress (Supplementary Table S7). To further examine differential rhythmicity within Type 2 ortholog groups, ortholog groups that also contained arrhythmic A. thaliana genes were filtered out, resulting in 2,258 ortholog groups that only contained rhythmic A. thaliana and arrhythmic K. fedtschenkoi genes. Clustering of the arrhythmic K. fedtschenkoi genes did not display a circular network typical of rhythmic genes, further validating the rhythmic gene expression detection method used in this study (Figure 3A). A slightly higher number of rhythmic A. thaliana genes were found with phase calls occurring during the night (Figure 3B). Within these 2,258 ortholog groups, 178 rhythmic A. thaliana TFs were identified and were phased to all phases of the day and night (Figure 3B).

Type 3: Rhythmic K. fedtschenkoi Genes With Arrhythmic Orthologs in A. thaliana
A total of 2,500 K. fedtschenkoi rhythmic genes were found to have no rhythmic A. thaliana orthologs in their respective ortholog groups (Figure 1C). Numerous processes were enriched in the K. fedtschenkoi gene set, including several metabolic processes and transcription from plastids (Supplementary Table S7). Differential rhythmicity was further examined and resulted in 1,759 ortholog groups that only contained rhythmic K. fedtschenkoi and arrhythmic A. thaliana genes. Arrhythmic A. thaliana genes did not cluster into a circular pattern as their rhythmic K. fedtschenkoi orthologs did ( Figure 3C). Plotting the distribution of phase calls of the rhythmic K. fedtschenkoi genes revealed a concentration of genes around midday ( Figure 3D). Within this set of ortholog groups, 233 K. fedtschenkoi TFs were identified, which were phased to all phases of the day and night and displayed a bimodal distribution of phase calls ( Figure 3D).
Type 4: Rhythmic A. thaliana Genes Only Without Orthologs in K. fedtschenkoi 788 ortholog groups were identified that had at least one rhythmic A. thaliana gene, totaling 2,293 rhythmic genes in A. thaliana, with no putative orthologs in K. fedtschenkoi (Figure 1C). Only toxin metabolic process was enriched in this gene set (Supplementary Table S7). There were 2,514 predicted rhythmic genes in K. fedtschenkoi with no detected orthologs in A. thaliana (Figure 1C). Only monoterpenoid biosynthetic process was enriched in this gene set (Supplementary Table S7). To determine if any genes in this gene set could be CAM-specific (i.e., shared by multiple CAM species but not by other non-CAM species) or K. fedtschenkoi-specific (i.e., not shared by other CAM and non-CAM species), protein sequences for each gene were BLASTed against a collection of plant proteomes spanning monocots and dicots and including all photosynthetic types (i.e., C 3 , C 4, and CAM) (Supplementary Table S3). Results indicate that 1,486 K. fedtschenkoi genes in this gene set were found in both non-CAM and CAM species (1 < e-10; Supplementary Table S10). Four K. fedtschenkoi genes were found to have homologs in only other CAM species and 443 K. fedtschenkoi genes were found to have homologs in only the non-CAM species. Of the four genes found in only CAM species, none were found in all of the six CAM species tested. Finally, 581 genes were found to be K. fedtschenkoi-specific. A majority of the rhythmic K. fedtschenkoi-specific genes had phase calls during midday, approximately 8 h after the beginning of light period (Supplementary Table S11), which coincided with the four CAM-specific K. fedtschenkoi genes (Supplementary Table S11).

Core Clock Genes in A. thaliana and K. fedtschenkoi
Investigating gene expression between A. thaliana genes and their orthologs in K. fedtschenkoi revealed cases of loss and gain of rhythmicity as well as phase shifts. One mechanism that can cause differences in gene expression between orthologs is differences between the regulatory network an ortholog is connected to. For example, through evolutionary processes, an ortholog can become connected to a different TF(s) in the network or become connected to an entirely different regulatory network, via reconnection with a new TF(s) or the its original TF(s) connected with a new network (Nadimpalli et al., 2015;Mateos et al., 2017). Both of these mechanisms can result in altered gene expression of an ortholog. The core circadian clock, which is a small regulatory network of interacting TFs, is known as a core mechanism that drives rhythmic gene expression in plants (Nohales and Kay, 2016). Therefore, we investigated the gene expression of core circadian clock TFs to determine if the loss and gain of rhythmicity and phase shifts seen above could be a result of changes in the core circadian clock network, e.g., phase shifts of or loss of rhythmicity in gene expression of core circadian clock TFs. Commonly, the A. thaliana core clock model (Nohales and Kay, 2016) is used to infer conservation of clock genes in other plant species. Following this, orthologs of the clock genes in A. thaliana were found in K. fedtschenkoi (Table 2 and   Supplementary Table S12). Variation in the K. fedtschenkoi core clock network was observed when assaying each component's expression dynamics. Core clock genes are typically characterized by cardinal circadian parameters, such as high amplitudes and fold-changes, along with highly statistically significant rhythmicity (Doherty and Kay, 2010;Zhang et al., 2014;Hughes et al., 2017). Using high amplitudes (max expr -min expr ) (FPKM units) as absolute amplitudes > 10, high fold-changes (max expr /min expr ) (FPKM units) as fold-changes > 2 and statistically rhythmic as JTK p ≤ 0.05, only 10 of the 23 K. fedtschenkoi core clock orthologs were found ( Table 2). Within these 10 genes, there were paralogs of RVE6 and ELF4. Interestingly, ELF3, ELF4, and LUX, which make up the evening complex (Nusinow et al., 2011) of the circadian clock network, displayed a concerted phase shift of 4 h ahead of their A. thaliana orthologs (Figure 4A). Of the remaining 13 core clock orthologs, four did not pass the gene expression threshold (max FPKM > 1) used before the analysis of gene expression rhythmicity ( Table 2).
Constructing the core clock model for K. fedtschenkoi with this information could suggest a different architecture of the clock ( Figure 4B).

DISCUSSION
CAM plants arose independently from diverse C 3 ancestors and developed two distinguishing features, strict temporal control of metabolism (e.g., CO 2 fixation, sugar accumulation) and temporally inverted stomatal movement patterns relative to C 3 plants. The temporal significance of these features has led to the  hypothesis that these features developed under the control of the circadian clock (Boxall et al., 2005;Cushman et al., 2008). To further test this hypothesis, we first identified rhythmic genes in a CAM and C 3 species and identified orthologs between the two species. We next examined the expression profiles looking for evidence of phase shifts between rhythmic orthologs in these two species. In this study, we revealed diversification of circadian rhythmicity between CAM and C 3 photosynthesis species. Moreover, differences were found between the K. fedtschenkoi core clock network and the A. thaliana model core clock network (Table 2), providing molecular evidence supporting the hypothesis that the core clock could have impacts on CAM evolution. Approximately 30% of all genes were found to be rhythmically expressed in both species, which is consistent with previous circadian studies grown under similar conditions and using whole-leaf samples (Michael et al., 2008;Filichkin et al., 2011). Prior studies have examined rhythmic transcriptomes of various plant species and among them, similar functional groups have been found enriched in their plant's respective rhythmic gene set (Harmer et al., 2000;Filichkin et al., 2011;Koda et al., 2017). Surprisingly, there was no overlap in enriched functional groups when comparing the entire A. thaliana and K. fedtschenkoi rhythmic gene sets to each other. Our analysis cannot explicitly explain why there is no overlap in enriched functional groups. However, a potential explanation could be that the circadian clock controls alternate biological processes in K. fedtschenkoi in response to extreme environmental fluctuations. For example, CAM plants typically experience more drastic temperature differences between night and day than their C 3 counterparts. Additionally, temperature perception by the clock has been found to be different between plants resulting in different regimes of genes being transcribed and therefore different biological processes to occur (Filichkin et al., 2011). It would be worth exploring how C 3 plants respond at the transcriptional level to similar temperature conditions to determine if any parallels exist which could give insight into the CAM-to-C 3 transition.
To determine if an increase in the temporal scale could further reveal similarities and differences in rhythmic biological processes, we examined phase specificity of biological processes. Photosynthesis and carbohydrate metabolism are plant processes that take place during the daylight hours and the circadian clock is known to play a crucial role in their regulation (Adams and Carre, 2011;Filichkin et al., 2011;Koda et al., 2017). Consistent with this, starch metabolism and photosynthesis were both enriched during midday for both species, in turn confirming conservation of these two pathways across C 3 and CAM species. Though the similarities cease there and the idea of the clock controlling alternative biological processes is further supported as several processes are phase-specifically enriched in K. fedtschenkoi and not in A. thaliana. For instance, organic acid synthesis was enriched in the rhythmic transcriptome of K. fedtschenkoi, in line with the strict schedule of CAM-related metabolism to prevent futile cycling between dark-period CO 2 fixation to malate and light-period malate decarboxylation (Borland et al., 2014). Furthermore, processes related to gluconeogenesis were found enriched in K. fedtschenkoi's rhythmic transcriptome and not in A. thaliana's rhythmic transcriptome. Wai et al. (2017) have recently shown that metabolic processes such as glycolysis and gluconeogenesis have day-night rhythms in the CAM plant pineapple. These results support the idea that maintaining metabolic homeostasis via strict scheduling of the associated genes is a defining feature of CAM plants.
Several studies have compared rhythmic gene sets between different species and found conserved core clock networks but divergent clock output networks (Keegan et al., 2007;Michael et al., 2008;Filichkin et al., 2011;Eckel-Mahan et al., 2013;Boyle et al., 2017;Yin et al., 2018). Our comparative analysis of rhythmic gene sets between A. thaliana and K. fedtschenkoi also revealed divergence in clock output networks, but also in the core clock network (Figure 4). Focusing on clock output first, two clock-related features, which were observed in Type 1 and 2 ortholog groups, account for the divergence in clock output. The first clock-related feature is of synchronization of rhythmic orthologs. These synchronized orthologs were enriched with functional groups related to photosynthetic processes and a large portion of these orthologs were found phased during the day, coinciding with when photosynthesis takes place in both species (Brautigam et al., 2017). These results are comparable with the results of phase-specific enrichment of functional groups across the entire rhythmic gene sets of both species, providing further support that the circadian clock plays an important role in photosynthesis in both C 3 and CAM species.
The second clock-related feature is conserved rhythmicity between orthologs but divergence in timing of expression. The TFs in Table 1 represent potential mechanisms by which the clock alters rhythmic gene expression in the context of CAM evolution beyond just CO 2 fixation. Specifically, the listed TFs could be associated with different core clock genes in K. fedtschenkoi than in A. thaliana, thus shifting not only their own gene expression but the genes they regulate. Several key CAM-related CO 2 fixation genes were found in the Type 1 ortholog groups that displayed this clock phenotype, although, none were found in the two groups of phase-shifted orthologs and TFs. However, this does not rule out the remaining TFs displaying a different phase shift potentially being the mechanism behind the phase shift of these key CAM-related CO 2 fixation genes. Regardless, the TFs in Table 1 warrant further functional investigation into how they are integrated into the output of the clock between species and what downstream processes they regulate.
The core clock network is one of the main causes of rhythmic gene expression seen in plants and has typically been found conserved across plant species (Filichkin et al., 2011;Koda et al., 2017;Sharma et al., 2017). In K. fedtschenkoi, several genes orthologous to A. thaliana core clock genes were found conserved in copy number and expression dynamics (Figure 4). However, some differences were observed between components of each species' core clock network. For example, the evening complex, consisting of ELF4, ELF3, and LUX, was found to be altered in gene copy number (i.e., ELF4) and in timing of gene expression (Figure 4). Variation in gene copy numbers between orthologous core clock genes is not surprising as it has been observed in other species, the mechanism likely being fractionation and/or loss of genes during speciation events (Lou et al., 2012). The shifting of expression in the components of the evening complex, mostly to midday (Figure 4), does present an inserting case, as a common theme seen in this comparative analysis was rhythmic K. fedtschenkoi genes, including TFs, being highly phased to midday (Figures 2, 3). Whether the shift in the evening complex is cause for the joint shift in rhythmic genes seen here cannot be determined in this study. However, the fact that other core clock genes did not display a similar shift in gene expression and rhythmic gene expression is typically driven by the core clock does give reason for further investigation into this potential relationship. Studies examining potential relationships between TFs in Table 1 and components of the evening complex would provide better insight into how the circadian clock played a role in the evolution of CAM.
For a further example, some components, such as CCA1 and TOC1, were found to have low expression ( Table 2), which is not typical of core clock genes (Lou et al., 2012;Yeung et al., 2018). This study is unable to identify the mechanism(s) behind the low expression values of these core clock genes; however, a possible cause could be feedback from clock output networks back into the core clock network. Feedback regulation between the clock and input signals as well as between output networks and the clock is a common feature of the circadian system (Fankhauser and Staiger, 2002;Li et al., 2011;Wenden et al., 2011;Chow et al., 2014;Kolmos et al., 2014). The low expression of some core clock genes could be a gating and/or compensation mechanism by the clock to some unknown factor(s). Alternatively, these plants were grown in diel conditions, so the impact of environmental variations, such as light, cannot be exclude as causes. To further investigate differential rhythmicity of K. fedtschenkoi core clock components with their respective A. thaliana orthologs, it will be necessary to generate time-course transcriptomesequencing data for leaf samples collected from both A. thaliana and K. fedtschenkoi at 2-h intervals during 48-h continuous light or dark period. A comprehensive analysis of time-course transcriptome-sequencing data, such as those described in Michael et al. (2008) and Filichkin et al. (2011), from A. thaliana and K. fedtschenkoi leaves under various light and temperature conditions would provide further insight in circadian control of CAM. Additionally, it was recently hypothesized that posttranscriptional regulation plays an important role in circadian clocks because a transcript does not oscillate does not mean that the protein levels are not rhythmic (Kojima et al., 2011;Lim and Allada, 2013). Future studies on rhythmic profiles of protein expression will be needed to gain a comprehensive understanding of circadian rhythm in plants.

DATA AVAILABILITY
All datasets generated for this study are included in the manuscript and the Supplementary Files.