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 (Yang et al., 2017). 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).

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 co-regulators 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).

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.

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).

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.

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 processes, photosynthesis (light harvesting), glyceraldehyde-3-phosphate 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.

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₂ fixation in K. fedtschenkoi (Yang et al., 2017) 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.

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).

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).

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₃, 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).

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).