High-Resolution Profiling of a Synchronized Diurnal Transcriptome from Chlamydomonas reinhardtii Reveals Continuous Cell and Metabolic Differentiation

The Majority of the Chlamydomonas Transcriptome Is Differentially Expressed in Synchronous Diurnal Growth Conditions

We optimized growth conditions and culture synchrony in order to identify genes that exhibit periodic diurnal and/or cell cycle-regulated expression patterns. Phototrophic cultures were grown in a 12-h-light/12-h-dark diurnal cycle under turbidostatic control to maintain uniform illumination, biomass density, and nutrient levels during the light phase, during which individual cells grew in mass by ∼10-fold. Immediately following the dark transition the cells divided synchronously on average three or four times to produce 8 or 16 daughters. Cultures were independently sampled over two 24-h periods at 1-h intervals beginning at ZT1 (Zeitgeber time = 1 h following the onset of illumination). Additional samples were taken at 30-min intervals between ZT11 and ZT15 during S/M phase, for a total of 28 samples per time course, ending at ZT24 (Figure 1A). Culture synchrony was assessed by comparing cell growth in replicate experiments (Figure 1B) and by monitoring progression through two cell cycle transitions: Commitment (reached after the minimum amount of growth necessary to complete one division cycle) and S/M (Figure 1C). The average cell density during the light phase decreased as the average cell size increased so that uniform biomass density and illumination were maintained throughout the growth phase (Figure 1B).

• Download figure
• Open in new tab
• Download powerpoint

Figure 1.

Experimental Design and Culture Metrics.

(A) Clock diagram showing experimental design. Diagram shows ZT time of light (ZT0 to ZT12) and dark (ZT12 to ZT24/0) periods of synchronized Chlamydomonas culture with growth (G1), division (S/M), and resting (G0) stages marked by inner gray arrows. Cartoons show relative cell sizes and two successive divisions by multiple fission. Cultures were sampled hourly with additional samples taken at 0.5-h intervals during cell division (brown lines).

(B) Cell size distributions from replicate cultures taken at 3-h intervals during the light phase.

(C) Plots showing fraction of each replicate culture that had passed Commitment and mitotic index at each time point.

(D) Heat map depicting correlation between replicates for the most variably expressed genes (coefficient of variation ≥ 1.2) during the diurnal cycle.

RNA was prepared from the samples and processed for high-throughput transcriptome sequencing to generate more than three million uniquely mapping reads per sample (Supplemental Table 1). The subset of genes that showed the most variable expression levels over each diurnal cycle (coefficient of variation ≥ 1.2) was used to calculate Pearson correlation coefficients of all pairwise combinations between replicate samples (Figure 1D). Replicates were highly correlated (R ≥ 0.97) and could be subdivided into three major groups based on intersample correlation values: (1) light phase (ZT2-ZT12), (2) light-dark transition (ZT12.5-ZT15, and (3) dark phase (ZT16-ZT1) (Figure 1D).

Of the 17,737 predicted genes in Chlamydomonas (Merchant et al., 2007; Blaby et al., 2014), 15,325 had read abundances above a minimal cutoff of 1.061 RPM (reads per million mapped), which was chosen as an inclusion threshold based on reliability for detecting differential expression (see Methods). Of these, 12,592 showed at least a 2-fold change from their mean value in one or more time points across the diurnal cycle with a false discovery threshold of 0.05 and a maximum abundance >1 RPKM (reads per kilobase per million mapped). These genes, representing >81% of the measurable transcriptome, were grouped into 18 major clusters (c1 to c18) based on similarity of expression patterns identified by the K-means algorithm (Soukas et al., 2000) implemented in the MeV software package (Saeed et al., 2003; see Methods), with each cluster containing 250 to 1071 genes (Figure 2A; Supplemental Data Set 1).

• Download figure
• Open in new tab
• Download powerpoint

Figure 2.

Coexpression and Functional Annotation of Genes with Diurnally Cycling Expression.

(A) Expression heat map of significantly differentially expressed genes. Expression of 12,592 genes during the diurnal cycle in 18 clusters (c1 to c18) is plotted as the number of standard deviations (σ) from the mean. Relative time, light and dark phases, and cell cycle stages are indicated above (Supplemental Data Set 1).

(B) MapMan term enrichment for clusters c1 to c18. Data are plotted as z-scores with high scores indicating significant enrichment. Level 1 and level 2 classifications separated by a comma are shown for selected terms from Supplemental Data Set 3.

We also identified genes with the least variable expression to find those that might be useful for normalization controls. A total of 171 genes had profiles with <2-fold expression variation between maximum and minimum RPKM (Supplemental Data Set 2) and were further examined in 15 other published RNA-Seq experiments (González-Ballester et al., 2010; Castruita et al., 2011; Fang et al., 2012; Urzica et al., 2012b; Malasarn et al., 2013; Goodenough et al., 2014; Schmollinger et al., 2014) (Mn- and Cd- experimental data; S.S. Merchant, unpublished data). We found 24 genes in our low-variance set with maximum fold expression changes of <2.5 and seven genes with fold expression changes of <2 in all other experiments. These stably expressed genes are predicted to be reliable internal controls for gene expression measurements under diverse conditions.

We compared our clustering method to another algorithm, JTK (Jonckheere-Terpstra-Kendall) cycle, which was designed for detecting rhythmic expression from transcriptome data (Hughes et al., 2010). JTK cycle identified 12,534 rhythmic genes with a stringent false discovery threshold (1E-10), 11,142 of which overlapped with our periodic expression data set (Supplemental Figure 1A). In addition, the phasing (i.e., LAG) assignments from JTK cycle were well correlated with specific diurnal expression clusters, meaning that, overall, both methods identified similar periodic expression patterns (Supplemental Figure 1B). However, we also noted a sizeable fraction of genes with either complex or transient expression profiles that were either missed by JTK cycle or given low-confidence scores despite their high reproducibility between replicates and high expression amplitudes (maximum expression/mean expression) (Supplemental Figure 1C). Conversely, the genes identified by JTK cycle that were not included in our differential diurnal set were almost entirely those that cycle reproducibly, but with low amplitude (Supplemental Figure 1D). We therefore used our original 18 coexpression clusters as a primary basis for further evaluation.

Functional Annotation of Expression Clusters

To detect functional specialization within the 18 expression clusters, we tested for enrichment of MapMan annotation/ontology terms (Thimm et al., 2004). Term enrichments found in different clusters included protein synthesis (c2, c4, c5, and c12) photosynthesis (c4, c6, c5, and c7), cell division (c10), DNA synthesis and chromatin structure (c1, c9, and c10), abiotic stress (c2 and c7), and cell motility (c11 to c13) (Figure 2B; Supplemental Data Set 3). Many of these ontologies were pursued further using manually curated sets of genes that provided more accurate and comprehensive bases for analysis of coexpression.

Elucidation of Coordinate Expression Patterns and Phase Relationships between Cell Cycle and Chloroplast Division Genes

We extended previous studies of cell cycle gene expression (Bisová et al., 2005; Fang et al., 2006) by first performing comprehensive annotation for 108 Chlamydomonas genes predicted to be involved in cell cycle regulation, DNA replication, chromosome segregation, and chloroplast division (Supplemental Data Set 4). Entry into S/M phase began just after ZT11 and was largely complete by ZT15 (Figure 1C). Strikingly, nearly every one of the cell cycle-related genes showed coordinated upregulation and high peak amplitude just prior to or during S/M phase, with a strong enrichment of replication genes in cluster c10 (Figures 2B and 3A to 3D; Supplemental Figures 4 and 5 and Supplemental Data Set 4). Although individual mother cells do not cycle through each round of S phase and mitosis in synchrony (Fang et al., 2006), phasing of cell cycle gene expression is informative because the first genes expressed will correspond to the earliest cell cycle events (e.g., S phase), while genes that are expressed later will correspond to subsequent events (e.g., mitosis or chloroplast division).

• Download figure
• Open in new tab
• Download powerpoint

Figure 3.

Coordinate Expression and Phasing of Cell Cycle Genes.

Plots of relative expression during light and dark phases. Phases are indicated by white and gray shading with expression levels normalized to a maximum of 1. Replica-averaged expression profiles of individual genes (italicized labels) or averages for all members of functional groups (nonitalicized labels) within each category are shown. Plots were smoothed as described in Methods. Additional expression data for cell cycle genes are in Supplemental Figures 4 to 11 and Supplemental Data Set 4.

(A) to (C) DNA replication-related genes and protein complexes.

(D) Structural maintenance of chromosomes (SMC) complexes.

(E) Predicted S phase and mitotic cyclins and CDKs.

(F) Additional cell cycle regulators.

(G) APC/C genes.

(H) RB tumor suppressor pathway genes.

(I) D-type cyclin genes.

(J) Chloroplast division genes.

Forty-six DNA replication-related genes showed coordinated expression, peaking at approximately ZT11 near the beginning of S/M, and this coordinated expression serves as a useful temporal landmark for other cell cycle events (Figures 3A to 3D; Supplemental Figures 4 to 6). Cyclin-dependent kinases CDKA1 and CDKB1 are homologs of universally conserved or green lineage-specific cell cycle regulators, respectively (Robbens et al., 2005; Bisová et al., 2005; De Veylder et al., 2007; Cross and Umen, 2015). In Chlamydomonas, CDKA1 function is required to initiate DNA replication, whereas CDKB1 function is required later for mitotic progression (Tulin and Cross, 2014). CDKA1 expression peaked at the time of DNA replication but its lowest expression during G1 phase never dropped below 20% of its peak value (Figure 3E), consistent with CDKA1 protein being detectable in G1 phase cells (Oldenhof et al., 2004). In contrast, CDKB1 expression is near zero when cells are not dividing and peaked sharply at ZT12, an hour after CDKA1 (Figures 3E; Supplemental Figure 7B). Two putative mitotic/S phase cyclin genes, CYCA1 and CYCB1, were expressed maximally at ZT13, just after CDKB1, while the hybrid A-B cyclin gene CYCAB1 was expressed maximally at ZT9, 2 h prior to visible signs of cell division (Figures 3E; Supplemental Figure 7B). These data suggest that CYCAB1 acts as an early activator of cell division and/or S phase, while CYCA1 and CYCB1 act later as mitotic cyclins. Other conserved genes with predicted roles during mitosis (WEE1, ESP1, and CKS1) were expressed in a pattern similar to CDKB1 (Figures 3F; Supplemental Figure 7C and Supplemental Data Set 4).

The anaphase promoting complex/cyclosome (APC/C) is an E3 ubiquitin ligase that triggers mitotic exit and also serves to maintain low CDK activity in G1 phase (Morgan, 2007; Chang and Barford, 2014). Genes encoding APC/C subunits were expressed in a broad peak from ZT10 to ZT13 but were also expressed detectably in G1 phase (Figures 3G; Supplemental Figure 9A) as was the gene encoding the G1 phase APC/C activator CDH1. In contrast, the gene for the predicted mitotic activator of APC/C, CDC20, was expressed in a sharp peak around ZT13 with little or no expression outside of S/M phase (Figures 3G; Supplemental Figure 9B). These data match key findings in other species where there is a highly specific mitotic requirement for APC^CDC20 and a broader requirement for APC^CDH1 upon mitotic exit and during G1 phase (Zachariae and Nasmyth, 1999).

Proteins of the retinoblastoma (RB)-related protein complex in Chlamydomonas are present throughout the cell cycle where they are required for cellular size control at Commitment and to regulate cell division number during S/M phase (Umen and Goodenough, 2001; Bisová et al., 2005; Fang et al., 2006; Olson et al., 2010). mRNAs for each subunit of the complex (MAT3 [RBR], E2F1, and DP1) are detectable in all samples with peak expression for all three occurring relatively early in the cell cycle at ZT9, consistent with a role in initiation of S phase and/or earlier events such as Commitment (Figures 3H; Supplemental Figure 7A). Interestingly, a gene encoding a E2F-related protein, E2FR1, which has degenerate DNA binding and dimerization domains (Bisová et al., 2005), is expressed prior to the genes encoding the core RB complex, suggesting that this protein has a function that is distinct from that of the canonical RB complex (Figures 3H; Supplemental Figure 7A).

Besides recapitulation of known or anticipated regulatory patterns, our data provide a means for classifying groups of cell cycle regulators that have been less extensively investigated in Chlamydomonas. Whereas most Chlamydomonas cell cycle genes are present in a single copy (Bisová et al., 2005), the D-cyclin family has four paralogs (encoded by CYCD1-CYCD4) whose origins are at least as old as the split between Chlamydomonas and its multicellular relative Volvox carteri (Prochnik et al., 2010). Each D cyclin has a distinct expression profile with peaks at ZT8 (CYCD4), ZT10-ZT11 (CYCD2), ZT12 (CYCD1), and ZT15 (CYCD3) (Figure 3I; Supplemental Figure 7D). These divergent expression patterns point toward D-cyclin subfunctionalization associated with their expression peak phasing: CYCD4 is a candidate regulator of and/or marker for Commitment; CYCD2 is a candidate S phase activator; CYCD1 and CYCD3 are candidate mitotic activators; and the prolonged expression of CYCD3 and CYCD2 also suggests a postmitotic function for these genes.

Unlike land plants that have many individual chloroplasts per cell, Chlamydomonas and most other unicellular Chlorophyte algae have a single large chloroplast that must be physically partitioned during cell division, a process that occurs before cytokinesis (Goodenough, 1970; Gaffal et al., 1995) (Figure 4). Several genes predicted to encode chloroplast division proteins in Chlamydomonas were previously shown to be expressed during S/M (Wang et al., 2003; Adams et al., 2008; Hu et al., 2008; Miyagishima et al., 2012) but were not sampled at high temporal resolution. We found that known chloroplast division genes peak at ZT11-ZT12, consistent with chloroplast division preceding mitosis and cytokinesis, whose cognate genes are expressed at ZT13 (Figures 3J and 4; Supplemental Figure 10). A number of uncharacterized genes with predicted chloroplast-targeting sequences (364 genes) and/or membership in the GreenCut group encoding conserved plant/algal proteins (55 genes) (Karpowicz et al., 2011) were coexpressed with the known chloroplast division genes in clusters c10 or c11 and are candidates for direct participation in chloroplast division or other chloroplast-related processes that are coordinated with the cell cycle (Supplemental Data Set 4).

• Download figure
• Open in new tab
• Download powerpoint

Figure 4.

Summary of Cell Cycle Gene Expression.

Graphical representation of expression phasing for genes or gene groups related to cell division and cell cycle control (ZT11 to ZT14) with darkest colored shading indicating peak expression time. The lower cartoon depicts key stages of division and the relative timing of events within a single division cycle. The backward arrow indicates optional iterations of S phase and mitosis. Additional expression data for cell cycle genes are in Supplemental Figures 4 to 11 and Supplemental Data Set 4.

c10 membership was significantly overrepresented (z-score = 26) within our curated group of cell cycle-related genes identified primarily through homology searches (Figure 2B; Supplemental Figure 11 and Supplemental Data Sets 3 and 4). We asked whether membership in c10 also correlates with cell cycle function in genes identified in an unbiased genetic screen for temperature-sensitive lethal cell cycle mutants (Tulin and Cross, 2014). We found that genes corresponding to the DIV class of mutants, whose defect was in mitotic progression, were highly enriched for membership in c10 (z-score = 13) (Supplemental Figures 11 and 12 and Supplemental Data Set 5). In contrast, genes for a second class of mutants called GEX that showed a failure to exit G1 phase were found in several clusters and only showed significant enrichment in c1 (z-score = 3.76), which contains a large number of stress-related and chromatin-related genes (Figure 2B; Supplemental Figures 11 and 12 and Supplemental Data Sets 3 and 5). We predict from these data that DIV genes will likely have functions directly associated with cell cycle progression (e.g., DNA replication, mitosis, and cytokinesis), whereas GEX genes represent a more pleiotropic set of functions than DIV genes, and their further study may reveal new connections between the cell cycle and stress pathways and/or chromatin dynamics.

Expression of Genes Encoding Flagella and Basal Body Proteins Is Coordinated with the Cell Division Cycle

The flagella of Chlamydomonas are motility and sensory organelles that are homologous in structure and function to animal cilia and to the cilia or flagella found in other eukaryotic taxa (Carvalho-Santos et al., 2011; Ostrowski et al., 2011). During the cell cycle, flagella/cilia are resorbed or severed prior to S phase or mitosis, thereby freeing basal bodies to act as centrioles during mitosis and cytokinesis (Johnson and Porter, 1968; Coss, 1974; Plotnikova et al., 2009; Parker et al., 2010; Kobayashi and Dynlacht, 2011). Flagella are reformed upon mitotic exit (Wood et al., 2012) and also play a role in daughter cell release as a secretory site for hatching enzyme (Kubo et al., 2009). In our experiment, flagella resorption occurred synchronously as cells entered S/M and flagella reformation was scored upon hatching, by which time daughters already had full length flagella (Supplemental Figure 13).

Although a selected set of genes relating to flagella function were found to be upregulated just after cell division genes (Wood et al., 2012), we wanted to conduct a comprehensive analysis of basal body and flagella gene expression during the cell cycle. We first curated a core set of 193 basal body and flagella genes (Basal Body, Axoneme, IFT, and BBsome) that have been identified in previous studies (Keller et al., 2005; Lechtreck et al., 2009; Dutcher, 2009; Shiratsuchi et al., 2011; Bower et al., 2013) (Figures 5A to 5C; Supplemental Data Set 6). Clusters c11 to c14 are significantly enriched for the MapMan ontology term “cell motility” (z-score > 2.5)(Figure 2B; Supplemental Data Set 3) and contain an overrepresentation of known flagella and basal body genes (P < 0.02) (Figures 6A; Supplemental Figure 14). Remarkably, nearly every flagella and basal body gene exhibited a similar expression pattern with an abrupt rise in expression around ZT12 when S/M was ongoing, peak expression around ZT13-ZT14, and tapering expression during the remainder of the dark period (ZT15 to ZT24) (Figures 5A to 5C; Supplemental Data Set 6). This expression pattern matches the expected requirement for peak flagella synthesis in postmitotic daughter cells. Also notable were the expression patterns of several flagella and basal body genes that differed from the majority (Figures 5A to 5C; Supplemental Data Set 6) and correlate with additional functions for the gene products not involving flagella, e.g., in the actin cytoskeleton (IDA5 encoding actin and PRO1 encoding profilin) and in stress responses (HSP70A, DNAJ1, and HSP90A) (Schroda and Vallon, 2009). Therefore, our data suggest that additional basal body/flagella genes that are not part of clusters c11 to c14 (e.g., POC13/FMO11, POC17/PHB1, and CCT3)(Figure 5; Supplemental Data Set 6) may have additional functions outside of these organelles and that coexpression is a useful discriminatory filter for identifying such genes.

• Download figure
• Open in new tab
• Download powerpoint

Figure 5.

Flagella and Basal Body Gene Expression.

(A) to (C) Heat maps depict relative expression levels for flagella and basal body genes that are divided into functional subgroups. The maximum expression level for each gene is set to 1. Cell cycle stages, diurnal cycle time, and light and dark phases are indicated above each heat map. Individual gene names are on the right. Cluster membership of genes is shown to the left (gray bars); brown stars indicate nondifferentially expressed genes.

(A) Flagella axoneme genes and subgroups.

(B) Core basal body genes, POC genes, and BUG genes (described in the main text).

(C) Intraflagellar transport (IFT) and Bardet-Biedl syndrome protein complex (BBSome) genes. IFT genes are further divided into Complex A subunits, Complex B subunits, anterograde motor subunits, and retrograde motor subunits.

(D) Diagram of flagella and basal body associated structures at anterior of cell with expanded cross section of an axoneme. Additional data on flagella and basal body genes are in Supplemental Figures 13 and 14 and in Supplemental Data Sets 6 and 7.

• Download figure
• Open in new tab
• Download powerpoint

Figure 6.

Flagella and Basal Body Gene Functional Classification.

(A) Cluster membership distribution of flagella gene groups. Fractional compositions of the expressed transcriptome and flagella- or basal body-related gene subsets including the known flagella and basal body genes corresponding to those in Figures 5A to 5C, CiliaCut genes, FAP genes, and deflagellation response genes. Enrichment for membership of each gene set in c11 to c14 compared with the entire transcriptome is indicated by asterisks showing significant P values.

(B) Venn diagrams show overlapping membership among indicated categories of flagella-related genes. Cilia Cut, FAP, deflagellation response genes belonging to either c11 to c14 (upper diagram) or all other clusters (lower diagram) were examined. The observed number of genes in a sector is shown inside the sector in bold with the expected membership number and sd based on a null model shown below each number in parentheses. Asterisks indicate sectors whose observed membership deviates significantly from the expected value. Additional data on flagella and basal body genes are in Supplemental Figures 13 to 15 and in Supplemental Data Sets 6 and 7.

Among the basal body genes (Figure 5B; Supplemental Data Set 6), we noted significant differences in expression patterns and cluster membership among three subgroups: core basal body (those with validated functions), POC (proteome of centriole), and BUG (basal body upregulated after deflagellation): The latter two groups were identified in a basal body proteomic study, with POC genes validated by homology with centriolar proteins of other species and BUG genes validated by upregulated expression after deflagellation (Keller et al., 2005). We observed that core basal body and axonemal gene groups had distinct distributions among c11 to c14 (P = 6.2 E-09, Fisher’s exact test), with core basal body genes predominantly found in c11 (peak expression at ZT12) and axonemal genes predominantly found in c12 or c13 (peak expression at ZT13 and ZT14) (Figures 5A and 5B; Supplemental Figures 14D and 14E). POC gene cluster membership was not detectably different from core basal body gene cluster membership (P = 0.76, Fisher’s exact test), whereas BUG gene cluster membership was different from both axonemal and core basal body genes, but more similar to the axonemal cluster distribution (P = 0.026, Fisher’s exact test) than to the basal body cluster distribution (P = 8.2 E-05, Fisher’s exact test) (Supplemental Figures 14D and 14E). These data are likely to reflect functional differences and temporal ordering between the core basal body and POC genes that are required during S/M for basal body assembly, replication, and mitosis, and the BUG genes that we predict will encode basal body proteins related to postmitotic nucleation and assembly of flagella (e.g., transition zone proteins). Indeed, many of the BUG genes are also present in the flagella proteome, whereas no POC genes are in this subset (Supplemental Data Set 7). Although the specific functions of most POC genes are unknown, the POC3/CEP290 gene product localizes mainly to the transition zone and is required for flagella assembly, and therefore does not match our prediction for involvement in basal body assembly or replication (Craige et al., 2010).

Three additional large data sets of candidate flagella/basal body protein coding genes were investigated based on their synchronous expression profiles: genes encoding proteins of the flagella proteome (FAPs) (Pazour et al., 2005), CiliaCut genes that are conserved in ciliated species but missing from species without cilia or flagella (Merchant et al., 2007), and an RNA-Seq-based transcriptome study of genes upregulated after deflagellation (Albee et al., 2013). We asked whether our coexpression data would discriminate between genes encoding proteins with core functions in flagella or basal bodies versus those with additional functions, or in the case of the deflagellation experiment, between genes whose expression is triggered by the stress of deflagellation and genes that are required for flagella biogenesis. First, we found that the uncharacterized members for all three large data sets described above (FAP, CiliaCut, and deflagellation upregulated) are significantly enriched for membership in c11 to c14 (Figure 6A; Supplemental Figure 14C) and predict that the c11 to c14 subset of genes from these data sets has functions directly related to flagella biogenesis.

We extended and validated this observation by examining comembership of genes in the three data sets (FAP, CiliaCut, and deflagellation upregulated) with the assumption that genes found in two or three of the groups are more likely tied to flagellum/basal body function than genes found in only one of the groups. Indeed, we found that members of c11 to c14 were highly overrepresented in the set of genes that belong to more than one data set (P = 1.5 E-27), while genes in other clusters were underrepresented in the overlap set (P = 3.0 E-27) (Figure 6B; Supplemental Figure 15B).

In summary, our data demonstrate robust, coordinated expression timing of basal body and flagella genes that is likely coupled to cell cycle progression and can be used as a powerful filter for identifying and classifying proteins that function in these organelles.

Ribosomal Protein Genes Are Expressed in Compartment-Specific Temporal Patterns

Green organisms must coordinate the activity of three separate protein translation systems located in the cytosol, chloroplast, and mitochondria respectively. rRNAs are encoded by the nucleus and respective organellar genomes, but ribosomal protein genes (RPGs) for all three compartments are nucleus-encoded. In principle, RPG expression for the three compartments could be globally coordinated and directly reflects patterns of diurnal cell growth, but instead we found that RPGs were expressed in distinct compartment-specific patterns (Figure 7; Supplemental Data Set 8). Cytosolic RPGs were expressed at high levels throughout the diurnal cycle used in this experiment and as a group constitute between 15 and 45% of all unique protein coding gene transcripts at any given time (Supplemental Figure 16). While cytosolic RPGs almost never dropped below half their peak expression values, they showed a consistent increase during the dark period with peak expression at around ZT15 when cells exited S/M and new daughters were hatching (Figures 7A and 7D). Notably, this was not a time when cells were actively growing. Mitochondrial RPGs showed an almost mirror-image pattern compared with cytosolic RPGs with highest expression in the light and reduced expression in the dark when respiration is most active (Figures 7C and 7D). Most strikingly, chloroplast RPGs showed a very strong periodic expression pattern, with a sharp peak early in the light period at ZT2 followed by a continuous decline through the remainder of the light period and partial recovery during the dark period (Figures 7B and 7D). These expression patterns suggest an unexpectedly strong partitioning of ribosome biogenesis among the three compartments, and they underscore the unanticipated coupling of ribosome biogenesis to organelle-specific requirements rather than cell growth as a whole.

• Download figure
• Open in new tab
• Download powerpoint

Figure 7.

Expression Patterns of Genes Encoding Chloroplast-, Mitochondria-, and Cytosol-Targeted Ribosomal Proteins.

(A) to (C) Heat maps depicting relative expression levels of RPGs. Relative expression levels of RPGs and translational regulators targeted to cytosol (A), chloroplast (B), and mitochondria (C). The maximum expression level for each gene is set to 1. Cell cycle stages, diurnal cycle time, and light and dark phases are indicated above each heat map. Gene names are on the right with known targets of chloroplast translational regulators in parentheses next to the corresponding gene name. Brown stars indicate genes without significant differential expression (<2-fold deviation from mean as defined in Methods).

(D) Graph of averaged compartment-specific RPG expression. Average absolute expression levels of large and small subunit RPGs during the diurnal cycle for cytosolic (blue), mitochondrial (orange), and chloroplastic (green) subunits. Additional data for ribosomal protein gene expression and translational regulator gene expression are in Supplemental Figure 16 and Supplemental Data Set 8.

Outside of the general patterns described above, we found that several RPGs that encode plastid-specific ribosomal proteins (i.e., PSRPs; Zerges and Hauser, 2009) are expressed out of phase with the majority of chloroplast RPGs. These include PSRP-1, whose expression peaks sharply and transiently at ZT13 but has very little expression at ZT2 during peak expression for most chloroplast RPGs. PSRP-1 is a homolog of the cyanobacterial gene lrt, whose transcript also accumulates after the light-to-dark transition. lrt is thought to stabilize ribosomal subunit association, possibly as a mechanism for stress-regulated translational control (Tan et al., 1994; Samartzidou and Widger, 1998), and the identification of the same pattern in Chlamydomonas suggests a deeply rooted common function for PSRP-1 in the green lineage. RAP38 and RAP41 are related proteins in the GreenCut group of conserved green lineage proteins (Karpowicz et al., 2011) that purify as stoichiometric components of intact Chlamydomonas 70S chloroplast ribosomes (Yamaguchi et al., 2003) and whose genes have nonoverlapping diurnal expression patterns (Figure 7B). Mutants of the Arabidopsis homologs of RAP38 and RAP41 (CSP41a and CSP41b, respectively), have defective rRNA processing and altered polysome profiles, indicating a role for this conserved pair of proteins in rRNA processing and ribosome biogenesis (Beligni and Mayfield, 2008). The noncanonical expression patterns we identified for RAP38 and RAP41 further suggest that they may function differently than core ribosomal proteins and possibly from each other.

In addition to RPGs, we examined expression patterns for other translational machinery and for nucleus-encoded regulators of chloroplast gene expression (Zerges and Hauser, 2009), some of which also showed atypical expression compared with chloroplast RPGs and suggest dynamic changes in chloroplast translation profiles at specific times of the diurnal and/or growth and division cycles (Figure 7B; Supplemental Data Set 8).

In summary, our data support a highly orchestrated pattern of ribosome assembly and compartment-specific translational control. Coexpression profiling further organizes translation-related genes into specific expression subgroups that will provide guidance for future studies of protein biosynthesis.

Plastid and Mitochondrial Protein Complex Subunits Show Phased Expression Patterns

Metabolism is driven by multisubunit protein complexes in chloroplasts and mitochondria, which generate ATP and reductant using either light energy (chloroplasts) or energy from catabolism (mitochondria). Light-driven reactions in the chloroplast are mediated by photosystems I and II (PSI and PSII) and associated light-harvesting complexes (LHCI and LHCII), along with cytochrome b₆f and the plastid ATP synthase complex. In mitochondria, the electron transport chain has five major complexes (I to V), including an ATP synthase complex that is encoded by a different set of genes than is the chloroplast ATP synthase. We examined gene expression patterns for the nucleus-encoded subunits for each of these complexes to determine both whether they are coordinately expressed and whether there are peak expression phase differences between them.

Subunits of each photosynthetic complex show broad light-phase expression peaks and are coexpressed in similar patterns, but they also display reproducible intercomplex phasing differences (Figures 8; Supplemental Figure 17 and Supplemental Data Set 9). The earliest expressed genes encode subunits of the ATP synthase and b₆f complexes whose mRNA abundances reach their peaks by ZT4-ZT6 and whose basal expression is relatively high compared with genes encoding photoactive complexes (Figure 8). Genes encoding subunits of PSI, PSII, and LHCI all show near zero expression in the dark and >1000-fold increases in transcript abundance in the light, with peaks at ZT6-ZT8, representing a 2-h phase delay compared with peak expression of genes encoding cytochrome b₆f and ATP synthase subunits. Genes encoding LHCII subunits showed a further phase delay with a peak from ZT8 to ZT10 (Figure 8; Supplemental Figure 17 and Supplemental Data Set 9). The staggered phasing we observed with cytochrome b₆f and ATP synthase complexes expressed earliest may enable productive photochemistry without delay upon PSI and PSII assembly, thereby minimizing phototoxic side reactions.

• Download figure
• Open in new tab
• Download powerpoint

Figure 8.

Coexpression and Phasing of Nucleus-Encoded Genes for Photosynthetic Complexes.

(A) to (F) Plots of relative expression of photosynthetic genes. Data are averaged for two replicates during light and dark phases (indicated by shading) with expression levels normalized to a maximum of 1 corresponding to the peak expression level of each gene.

(A) LHCII.

(B) LHCI.

(C) PSII.

(D) PSI.

(E) b₆f complex.

(F) CF₁-CF_o. Complete protein names encoded by each gene are in Supplemental Data Set 9.

(G) Composite expression patterns for each complex shown in (A) to (F). Data are smoothened as indicated in the Methods.

(H) Schematic showing the subunit architecture of each complex with color corresponding to graph colors gene in (A) to (F). Chloroplast-encoded genes are shown in pale gray. Additional data for photosynthetic complexes and assembly factors are in Supplemental Data Set 9.

Consistent with previous estimates of protein abundance (Oey et al., 2013; Drop et al., 2014; Natali and Croce, 2015), individual LHC gene mRNA abundances varied by several orders of magnitude compared with the peak expression levels of genes for PSI, PSII, b₆f, and ATP synthase complexes, which generally remained within a twofold range between genes (Figure 8; Supplemental Figure 17 and Supplemental Data Set 9). We also observed that genes encoding assembly factors for each photosynthetic complex were either coexpressed with or expressed just before the genes for their target complexes (Supplemental Figure 17 and Supplemental Data Set 9).

When grown photoautotrophically in a diurnal cycle, Chlamydomonas cells accumulate starch during the day through photosynthesis and catabolize it at night via glycolysis and respiration (Gfeller and Gibbs, 1984; Klein, 1987; Thyssen et al., 2001; Ral et al., 2006). The dark-phase respiratory activity associated with starch catabolism is reflected in the expression patterns of genes encoding subunits of mitochondrial electron transport/oxidative phosphorylation complexes I to V. Transcripts for these proteins began accumulating in the mid- to late-light period and peaked in the dark phase (ZT14 to ZT17) (Supplemental Figures 18 to 20 and Supplemental Data Set 10). Thus, respiratory complex genes are expressed out of phase with those of photosynthetic complexes and peak during the dark when respiration is expected to be most active.

Coordination of Tetrapyrrole Pathway Gene Expression with Light and Dark Phases

Chlorophyll and heme are the two major classes of tetrapyrroles in photosynthetic cells and both play key roles in energy metabolism. Genes for the common and chlorophyll-specific tetrapyrrole biosynthetic enzymes are coexpressed, with transcript levels increasing rapidly in the light and peaking at ∼ZT6. In contrast, HEM15, which encodes the ferrochelatase enzyme dedicated to heme biosynthesis, shows increasing expression during the dark, highest expression at the end of the dark period, and steadily decreasing expression during the light period when the chlorophyll pathway is most highly expressed (Supplemental Figure 21 and Supplemental Data Set 11). This temporal partitioning between heme and chlorophyll pathway genes may have evolved as a strategy to reduce competition between the two tetrapyrrole biosynthetic pathways and to ensure that chlorophyll is the major tetrapyrrole produced during the light phase (Supplemental Figure 21 and Supplemental Data Set 11).

Transient Stress Responses Occur at the Dark-to-Light Transition

The abrupt dark-to-light transition in our experiment enabled us to identify and refine a cluster of 280 genes, mostly derived from c1 and c2, whose expression is transiently induced in a sharp spike at the first light time point, ZT1 (Figure 9A; Supplemental Data Set 12) (see Methods). Functional annotation of this cluster showed enrichment for abiotic stress-related genes (Figure 9B) that include five heat shock protein (HSP) genes and a gene whose Arabidopsis ortholog (At1g04130) encodes a HSP cochaperone (Cre08.g375650) (Schroda and Vallon, 2009). The cluster also included VTC2, which encodes GDP-l-galactose phosphorylase, the enzyme that performs the first committed step in the Smirnoff-Wheeler pathway for ascorbate biosynthesis (Urzica et al., 2012a) (Figure 9A). We hypothesize that many other genes in this cluster are induced by the stress associated with the abrupt transition from darkness into high-light growth conditions.

• Download figure
• Open in new tab
• Download powerpoint

Figure 9.

Transient Stress Response at Dark-Light Transition.

(A) Transcript abundance for known stress response or ROS response genes (red) and GreenCut genes (green). The mean of all 290 light stress gene expression profiles comprising the light stress cluster is plotted in black. Light-gray shading indicates the dark period. Complete protein names encoded by each gene are in Supplemental Data Set 12.

(B) Functional characterization of the light stress cluster. MapMan ontology distributions with percentages of the total are shown in parentheses. Terms highlighted with a star are significantly enriched (P value < 0.05). In total, the light stress cluster includes 13 GreenCut genes, several of which have unknown functions (CGL150, CGL122, CPLD50, and CGL99) (Karpowicz et al., 2011; Heinnickel and Grossman, 2013). Additional data on the light stress cluster are in Supplemental Data Set 12.

LHC-like genes in the ELIP, LHCSR, and PSBS families are also associated with light stress responses (Heddad and Adamska, 2002; Hutin et al., 2003; Elrad and Grossman, 2004; Bonente et al., 2008). ELIP1-ELIP5 (encoding homologs of early light inducible proteins) have previously been described in Chlamydomonas (Elrad and Grossman, 2004; Teramoto et al., 2004), and, based on recent genome updates (Blaby et al., 2014), we identified five additional Chlamydomonas ELIP genes, ELIP6 to ELIP10 (Supplemental Figure 22 and Supplemental Data Set 13). ELIP2, ELIP3, ELIP4, ELIP9, and ELIP10 as well as one of two PSBS paralogs, the GreenCut member PSBS2 (Cre01.g016750), are in the light stress cluster. The other PSBS paralog, PSBS1 (Cre01.g016600), showed a similar expression pattern as PSBS2 (Supplemental Data Set 1), but its expression maxima was below the threshold used for inclusion in the cluster. The highly transient light-stress-activated expression of PSBS genes discovered in our study may explain their lack of reported expression in previous experiments and also suggests that, contrary to previous reports (Anwaruzzaman et al., 2004; Bonente et al., 2008; Peers et al., 2009), the Chlamydomonas PSBS proteins may play a role in nonphotochemical quenching (NPQ), the conversion and dissipation of excess excitation energy as heat.

Besides abiotic stress, other functional terms that are enriched in the stress cluster include ABC transporters, protein and amino acid biosynthetic genes, and protein translation (Figure 9B; Supplemental Data Set 14). This suggests that a concerted metabolic rewiring might occur at the dark-to-light transition when cells must rapidly respond to light stress, while also preparing for a massive upregulation of their photosynthetic capacity in the subsequent light period.

Central Carbon Metabolism Genes Show Pathway-Specific Expression Patterns

Central carbon metabolism in photosynthetic cells must respond to changes in growth conditions such as activity of photosynthetic complexes and external nutrients, and these changes can be reflected by coordinated changes in transcriptional networks for metabolic genes (Wei et al., 2006). The transcript abundance of genes encoding enzymes of the Calvin-Benson-Bassham (CBB) cycle of CO₂ fixation exhibited highly correlated expression patterns along with CP12, encoding a chaperone for GAP3, and the gene for RUBISCO ACTIVASE (Figure 10; Supplemental Data Set 15). Peak expression for nearly all genes encoding enzymes of the CBB cycle occurred between ZT5 and ZT8. An exception from this pattern was RBCS2, one of the two Rubisco small subunit paralogs that was expressed continuously in our diurnal growth conditions. Another exception was TPIC1, the single gene in Chlamydomonas encoding triose phosphate isomerase, which is expressed during the light and dark periods, possibly reflecting its dual roles in the CBB cycle as well as in glycolysis during the dark phase (Supplemental Data Set 15).

• Download figure
• Open in new tab
• Download powerpoint

Figure 10.

Expression of Central Carbon Metabolism Genes and Their Diurnal Phasing.

Metabolic pathways are shown for the CBB cycle (A), acetate metabolism (B), and the TCA cycle (C). Beside each pathway schematic are average expression estimates of the two replicates normalized to peak expression and raw RPKM values for each of the two replicates for each gene of the pathway. Gray shading indicates the dark period. The enzymes catalyzing each step are indicated with primary gene names beside each pathway arrow. Color coding (according to the heat map) and clock icons indicate the time of peak expression. Black arrows and enzyme names in black are used for steps where different isoforms or enzyme subunits are not expressed in the same pattern. Additional data on central carbon metabolism genes including protein names are in Supplemental Data Set 15.

Because our cultures were grown photoautotrophically (i.e., without acetate or other organic carbon sources), the patterns we observed for dark metabolism gene expression reflect a state where organic carbon is derived solely from endogenous pools such as starch or storage lipids. Genes encoding enzymes required for the tricarboxylic acid (TCA) cycle and acetate assimilation had less coherent expression patterns than those for the CBB cycle but were still strongly biased toward peak expression during the dark period (ZT12 to ZT24) when respiration dominates energy metabolism (Figure 10; Supplemental Data Set 15). Transcripts for starch biosynthesis enzymes gradually increased during the light period, whereas the expression of genes encoding starch catabolic enzymes was restricted to the dark period when starch is broken down (Supplemental Data Set 15).

Paralog Expression Clustering Predicts Pathway Assignments for Duplicated Metabolic Enzymes

Several enzymes of central carbon metabolism in Chlamydomonas are encoded as multigene families in which individual isoforms may be targeted to different subcellular compartments and are associated with pathways operating at different times during the diurnal cycle (Supplemental Data Set 15). Experimental localization of biochemical activities has been described for some enzymes (reviewed in Johnson and Alric, 2013), but the pathways within which many enzyme isoforms act in Chlamydomonas remain unknown. This uncertainty impacts the accuracy of metabolic models that try to simulate metabolism in a compartmentalized cell (Dal’Molin et al., 2011).

We used expression profiles of “signature” genes that are dedicated to specific pathways (e.g., sedoheptulose-1,7-bisphosphatase [SBP1] to the CBB cycle and isocitrate lyase [ICL1] to the glyoxylate cycle) as a template for correlation measurements that helped assign individual paralogs of duplicated genes to their cognate pathways (Supplemental Data Set 15). For example, citrate synthase activity is required for both the TCA cycle in mitochondria and for the glyoxylate cycle, whose location is presumably in microbodies (the algal equivalent of a plant peroxisome) (Hayashi et al., 2015). The expression profiles of two citrate synthase paralogs encoded by CIS1 and CIS2 were compared with those of the signature TCA cycle genes (SDH3 and SDH4) and to signature glyoxylate cycle genes (MAS1 and ICL1) using Pearson’s correlation. CIS2 expression was strongly correlated with the glyoxylate cycle genes (0.91) and more weakly with TCA genes (0.54), whereas CIS1 expression was poorly correlated with glyoxylate genes (−0.33) and moderately correlated with TCA genes (0.60) (Supplemental Data Set 15). The expression profile-based assignment of CIS2 to the glyoxylate cycle and CIS1 to the TCA cycle was supported and validated by data from Arabidopsis, in which the most similar homologs of CIS2 (CSY4 and CSY5) form a clade of peroxisome-targeted citrate synthases and the most similar homologs of CIS1 (CSY1 to CSY3) form a clade of mitochondrial-targeted citrate synthases (Pracharoenwattana et al., 2005). By applying a similar approach and an expression correlation threshold of ≥0.6, we propose pathway assignments for a total of 23 enzyme isoforms (Supplemental Data Set 15).

Transcription Factor Gene Expression Patterns

We examined the expression profiles of 230 predicted and/or functionally characterized DNA binding transcription factors (TFs) (Pérez-Rodríguez et al., 2010; Jin et al., 2014) in our study and from previously published transcriptome experiments (González-Ballester et al., 2010; Castruita et al., 2011; Fang et al., 2012; Urzica et al., 2012b; Malasarn et al., 2013; Goodenough et al., 2014; Schmollinger et al., 2014) (Mn- and Cd- experimental data; S.S. Merchant, unpublished data; Supplemental Data Set 16). We compared levels of induction, expression amplitudes, and absolute expression maxima as a way of categorizing the regulatory patterns of Chlamydomonas TFs. The mRNA abundances of 219 TFs were ≥1 RPKM in at least one time point of the time course, and of these, 191 showed significant differential expression patterns (Supplemental Figure 23 and Supplemental Data Set 16). Eighty-six differentially expressed TFs were not significantly upregulated in other experiments/conditions and are therefore possible candidates for controlling diurnal and/or cell cycle related gene expression. In contrast, 105 TFs showed differential expression in our experiment and significant upregulation in at least one other experiment, with 28 of these showing greater upregulation in other data sets compared with their peak amplitude in our data set (max_other/max_diurnal > 2) (Supplemental Data Set 16). Ten TFs were upregulated following deflagellation (Albee et al., 2013), with three of these also members of the core flagella clusters (c11 to c14) and, therefore, candidates for controlling expression of flagella-associated genes (Supplemental Figure 23 and Supplemental Data Set 16). In contrast, the other seven TFs upregulated in response to deflagellation were also induced in nutrient-limiting conditions, suggesting that they may be part of a more general stress response that occurs after deflagellation.

Poorly Expressed Genes and Condition-Specific Expression

The expression estimates for ∼2900 predicted genes in our data did not rise above 1 RPKM at any time point (Supplemental Data Set 17). It is likely that the products of these genes are required either at low levels or only under specific conditions that are not part of our diurnal regime (e.g., sexual reproduction and nutrient limitation). Alternatively, some among the 2900 poorly expressed genes may be pseudogenes. We compared expression estimates for these 2900 genes to those from other published Chlamydomonas transcriptomes (Supplemental Data Set 17) and applied a minimum expression estimate of 10 RPKM/FPKM as a threshold for condition-specific expression. In doing so, we identified a total of 460 genes that were significantly expressed in other transcriptomes and are therefore strong candidates for being condition-specific genes whose expression is not required for vegetative proliferation (Supplemental Data Set 17).