AraLeTA: An Arabidopsis leaf expression atlas across diurnal and developmental scales

Abstract

Mature plant leaves are a composite of distinct cell types, including epidermal, mesophyll, and vascular cells. Notably, the proportion of these cells and the relative transcript concentrations within different cell types may change over time. While gene expression data at a single-cell level can provide cell-type-specific expression values, it is often too expensive to obtain these data for high-resolution time series. Although bulk RNA-seq can be performed in a high-resolution time series, RNA-seq using whole leaves measures average gene expression values across all cell types in each sample. In this study, we combined single-cell RNA-seq data with time-series data from whole leaves to assemble an atlas of cell-type-specific changes in gene expression over time for Arabidopsis (Arabidopsis thaliana). We inferred how the relative transcript concentrations of different cell types vary across diurnal and developmental timescales. Importantly, this analysis revealed 3 subgroups of mesophyll cells with distinct temporal profiles of expression. Finally, we developed tissue-specific gene networks that form a community resource: an Arabidopsis Leaf Time-dependent Atlas (AraLeTa). This allows users to extract gene networks that are confirmed by transcription factor-binding data and specific to certain cell types at certain times of day and at certain developmental stages. AraLeTa is available at https://regulatorynet.shinyapps.io/araleta/.

Figure 1.

CIBERSORTx predictions for simulated bulk RNA-seq samples. A) Seventy-five percent of cells per cluster were used to train a CIBERSORTx model, while the remaining 25% of cells were repeatedly subsampled to generate 500 simulated bulk RNA-seq samples with mixed cell types (see the Materials and methods section). B) Pearson's correlation between the true cell-type proportions in the simulated bulk RNA-seq sample and the tissue proportions predicted by CIBERSORTx. All the correlations were statistically significant (P < 0.05), with Pearson's correlations >0.8 highlighted by the horizontal bar. C) The slope of the line of best fit for a comparison between the predicted and the true cell-type proportions (scatterplots in Supplementary Fig. S2). Values close to 1 (horizontal line) are most accurate. Taken together, these results suggest that CIBERSORTx can predict the relative proportions within the same cell type between RNA-seq samples, but that it consistently overpredicts or underpredicts certain cell types within an RNA-seq sample. D) Predicted cell-type proportions in microdissected leaf samples, with each row representing a leaf and each column representing a different predicted cell type. Note that phloem is short for phloem parenchyma.

Figure 2.

Shifts in cell-type proportions over developmental and diurnal time series. We predicted the proportion of cell types in a A) developmental time series (Woo et al. 2016) and B) diurnal time series (Hickman et al. 2017), utilizing a reference leaf scRNA-seq dataset (Procko et al. 2022). In A), each row represents a different RNA-Seq sample, with the proportions normalized using z-scores. The developmental stage and the biological replicate are indicated by the colored bars. In B), the proportions of cell types were predicted independently for each of the 4 biological replicates, and then these were averaged for each row and normalized by z-scores. The times are relative to dawn.

Figure 3.

The expression of light-sensitive genes in single-cell leaf data. Sets of genes were selected that are light-induced A) and B) or light-repressed C) and D) during either a nocturnal light treatment A) and C) or a light treatment after an extended night B) and D), based on Rugnone et al. (2013). The mean expression of these genes in cells in each cluster in Procko et al. (2022) was calculated, and z-scores were calculated for each column.

Figure 4.

Cell activity changes during bolting. A) We predicted the proportion of cell types in plants immediately before and after bolting (Redmond et al. 2023), utilizing a reference leaf scRNA-seq dataset (Procko et al. 2022). Each row represents a different RNA-seq sample (representing a single plant), with the proportions normalized using z-scores. B) The plants in A) are ordered by pseudotime, which was calculated by Redmond et al. (2023) as an arrangement of 65 individual plants sampled along a developmental trajectory. C) For each cell type, the Spearman ranked correlation was calculated between the proportion of that cell type in each plant versus another trait of the plant (biomass, leaf size, or pseudotime).

Figure 5.

Cell-type-specific gene expression of bolting-related genes. Using the high-resolution imputing function in CIBERSORTx, we predicted cell-type-specific expression of genes that were differentially expressed in bolted/unbolted plants. In this study, we show the results for the genes that have higher expression levels in bolted plants, with the inverse gene set shown in Supplementary Fig. S9. A) The z-score of the expression of these genes in the bulk RNA-seq experiment (Redmond et al. 2023), grouped by the cell type in which they were predicted to be expressed and ordered by pseudotime. B) The cell-type assignment, with red indicating that a gene is expected to be found in that cell type. C) The z-score of the mean expression of these genes in the scRNA-seq dataset (Procko et al. 2022).

Figure 6.

PAFway networks across different tissue types at different ages. Using PAFway, we generated networks of functional terms to assess changes in function between tissue types as they age. A) Mesophyll (Groups 1, 2, and 3) and vascular-associated (bundle sheath, phloem, vascular, companion, guard, and sieve) cells showed a large overlap in functional edges across age groups, with some network edges that are unique to specific tissues and ages. B) A network of functional terms that are associated with young mesophyll and vascular cells. C) A network of functional terms associated with old mesophyll and vascular cells. Together, these show the changing network over time, suggesting that these cells perform different physiological functions as they age. The colors of the edges in B) and C) correspond to the associated colored sections of the Venn diagram in A).

Figure 7.

A graphical abstract showing the cell-type-specific expression pattern. This refers to the relative changes in gene expression values across diurnal A) and developmental B) time series. The colors correspond to the color scale of the heatmaps in Fig. 2.