Age attenuates the transcriptional changes that occur with sleep in the medial prefrontal cortex

Abstract Sleep abnormalities are common with aging. Studies show that sleep plays important roles in brain functions, and loss of sleep is associated with increased risks for neurological diseases. Here, we used RNA sequencing to explore effects of age on transcriptome changes between sleep and sleep deprivation (SD) in medial prefrontal cortex and found that transcriptional changes with sleep are attenuated in old. In particular, old mice showed a 30% reduction in the number of genes significantly altered between sleep/wake and, in general, had smaller magnitudes of changes in differentially expressed genes compared to young mice. Gene ontology analysis revealed differential age effects on certain pathways. Compared to young mice, many of the wake‐active functions were similarly induced by SD in old mice, whereas many of the sleep‐active pathways were attenuated in old mice. We found similar magnitude of changes in synaptic homeostasis genes (Fos, Arc, and Bdnf) induced by SD, suggesting intact synaptic upscaling on the transcript level during extended wakefulness with aging. However, sleep‐activated processes, such as DNA repair, synaptogenesis, and axon guidance, were sensitive to the effect of aging. Old mice expressed elevated levels of immune response genes when compared to young mice, and enrichment analysis using cell‐type‐specific markers indicated upregulation of microglia and oligodendrocyte genes in old mice. Moreover, gene sets of the two cell types showed age‐specific sleep/wake regulation. Ultimately, this study enhances understanding of the transcriptional changes with sleep and aging, providing potential molecular targets for future studies of age‐related sleep abnormalities and neurological disorders.


RNA-Seq library preparation was performed by the University of Pennsylvania Next-Generation
Sequencing Core using a TruSeq Stranded mRNA kit (Illumina) and sequencing was done with 100 base-pair single-end reads on a HiSeq 4000 System (Illumina). On average (standard deviation), 24.5 (2.6) million total raw reads were obtained for each sample. Raw RNA-Seq reads were aligned to the mouse genome build mm9 by STAR version 2.5.3a [12]. After the alignment, an average of 17.3 (1.9) million reads were mapped uniquely to annotated genes. The aligned data were quantified at the gene level using scripts from the PORT pipeline (github.com/itmat/Normalization -v0. 8.4-beta). Fragments from the same strand as the annotated genes were identified as sense, and fragments from the opposite strand of the annotated genes were identified as antisense. Low-expressing genes were removed by keeping only the genes with mean log 2 counts per million (CPM) >-1.4 (equivalent to 5 counts with a library size of 15 million reads). A total of 14,656 genes with unique ensemble IDs that passed the filtering criteria were normalized using the "Trimmed Mean of M-values" (TMM) method in the edgeR package [13], and differential expression analyses were performed using the LimmaVoom package [14] (see also below). A multidimensional scaling (MDS) plot of all samples was generated using the plotMDS function in Limma with default parameters. The Euclidean distance between each pair of samples, calculated as the root-mean-square deviation of the top 500 most variable genes automatically selected for each pair of samples, was used as an indication of the similarity between samples [14]. The first dimension represents the leading-fold-change that best separates samples and explains the largest proportion of variation in the data, with subsequent dimensions having a smaller effect and being orthogonal to the ones before it. [15]

Assessment of Differential Expression
Gene expressions were measured at baseline (ZT0, 7AM, lights-on) in each age group (young, old) and at 4 time-points [ZT3 (10AM), ZT6 (1PM), ZT9 (4PM), ZT12 (7PM)] in each agegroup and behavioral state (SS, SD). Thus, there were a total of 18 experimental conditions defined by age and behavioral state. Using the LimmaVoom package in R software [14], we evaluated the differentially expressed genes (DEGs) between sleep deprived and spontaneous sleep behavioral states (SD vs. SS) within each age group. Differentially expressed genes between young and old (old vs. young) were only evaluated in the undisturbed conditions (SS and ZT0). For each of the 14,656 genes that passed through the filtering criteria, statistical significance evaluating differentially expressed genes between behavioral state or age-group was based on a moderated F-test [14,16] testing for any significant differences at any of the four time-points (e.g., an ANOVA-like test). Similarly, comparisons between young and old mice at baseline (ZT0) were performed using a moderated T-test (given a single time-point for comparison). Multiple comparisons correction for statistical significance was applied using the method of Benjamini and Hochberg [17] to control the overall false discovery rate (FDR) at 1%, unless otherwise noted.
Statistical interaction was tested between behavioral state and age group for the identified sleep genes and wake genes from both young and old animals combined. Analyses were performed in a similar manner to those described above, using a moderated F-test to assess whether the differences in gene expression between SS and SD differed between young and old mice at any of the four timepoints. Multiple comparison correction using Benjamini and Hochberg was used and genes with FDR less than 5% were considered significant for the interaction test; a more liberal FDR threshold was chosen given the expected reduced power to detect interaction effects.
Unless otherwise noted, the statistical tests comparing FDR values and fold changes of genes between young and old animals were performed using paired Wilcoxon signed-rank tests. Effect sizes were calculated using Cohan's d method [18].

Functional analysis
DAVID Bioinformatics Resources 6.8 was used for gene ontology (GO) enrichment analyses [19,20]. Unless otherwise noted, enrichment of functional annotation was performed against four databases: GOTERM_BP_Direct (biological process), GOTERM_MF_Direct (molecular function), UP_keywords (UniProt keyword), and KEGG pathways. The 14,656 genes used for generating DEGs were used as the background in all cases. Enrichment of individual terms was defined as a P<0.05 and containing at least 3 genes. Functional clustering of overlapping terms was performed with kappa similarity threshold >0.3, minimum similarity terms≥3, multiple linkage threshold >0.3, and EASE <0.05. The significance of the enriched functional clusters was evaluated using an enrichment score, which is equivalent to the -log 10 transformed geometric mean of the Ps of included functional annotation terms. Since only terms with P<0.05 are included in the clusters, the minimum enrichment score from enriched functional clusters is 1.3.

Gene set enrichment analysis
Competitive gene set enrichment tests were performed using the CAMERA [21] method in the Limma package [14]. The default setting of inter.gene.corr=0.01 was used for all tests.
Enrichment of the cell-type specific genes was done by matching customized gene sets containing lists of cell-type specific marker genes against the ranked fold-changes of different comparisons (e.g., SD vs. SS and Old vs. Young), as described above. Cell-type specific marker gene sets were obtained from three studies, a microarray study performed using isolated astrocytes, neurons, and oligodendrocytes from mouse forebrain by Cahoy et al. [22], a singlecell RNA-seq study performed using cells dissociated from mouse primary visual cortex by Tasic et al. [23], and a single-cell RNA-seq study performed with the use of a cocktail of inhibitors blocking neuronal activation caused by enzyme digestion during the dissociation process by Hrvatin et al. [24]. Table S9 lists the sources and selection criteria of the cell-type specific gene sets, as well as the full lists of the genes included in the gene sets. Figure S1. Venn diagram of the DEGs identified between SD vs. SS within the young and the old mice across timepoints (FDR<1%) after applying varying fold-change (FC) cutoffs (based on absolute log2FC averaged across four time-points). Numbers of DEGs passing FC-cutoffs decrease as the FC-cutoff threshold increases from zero (without FC-cutoff) to 0.1, 0.2, and 0.5 for both young and old mice. But the proportions of genes being common or age-specific between the young and old mice stay relatively unchanged. Figure S2. Variability of the gene expressions (counts per million) for the sleep and wake genes among biological replicates (n=6) was measured using coefficient of variance (CV) during sleep (left) and SD (right) and compared between young and old animals. (A) CVs are highly correlated between young and old animals at all the conditions tested. (B) stacked bar plots of the percentage of genes have increased CV in young (top) or old (bottom) mice. At only two of the eight condition (3h and 6h of SS) more than half of the genes (around 65%) showed increased CV in the old animals compared to the young animals. In the remaining six conditions (9h and 12h of SS, and SD at all four time points), 60% to 75% of the genes showed increased CV in the young animals than in the old animals. Figure S3. Scatter plots of the differences in -log 10 P-values of individual sleep (left) and wake (right) genes between old and young mice. Positive differences indicate increased significance in old mice, and negative differences indicate decreased significance in old mice. The insets show the bar plots of percentage of genes in each gene class, common or age-specific, that have decreased significance in old mice (negative differences) compared to the young mice. Specifically, more than 98% of the youngspecific sleep and wake genes, as well as 76.8% of common sleep genes and 67.2% of common wake genes are less significant in old mice, whereas 98.5% of the old-specific sleep genes and 94% of oldspecific wake genes are more significant in old mice. The small percentages of age-specific genes showing increased significance in the opposite age group are likely due to the applied fold-change cutoff that allocates similarly significant genes into age-specific groups. Combined togther, 76.2% of sleep genes (1624 of 2130) and 74% of wake genes (1589 of 2147) are less significantly regulated between sleep and wake in old animals compared to the young animals. Figure S4. Correlation between the delta -log 10 P-values and delta log 2 FC of individual sleep (red) and wake (blue) genes between old and young mice. Positive differences indicate increased significance or FC in old mice, and negative differences indicate decreased significance or FC in old mice. Combined pearson correlation coefficientce and p-value are shown. Figure S5. Percentages of expression changes relative to baseline (ZT0) of each gene during SS or SD were calculated and trend lines were plotted using Loess local regression (with 95% CI). Sleep genes were shown on the left (red headings) and wake genes were shown on the right (blue headings). Both common and young-specific genes showed a reduction in the magnitude of change between SS and ZT0 in the old mice compared to the young. Y.SS=young mice during SS; Y.SD=young mice during SD; O.SS=old mice during SS; O.SD=old mice during SD. Figure S7. Venn diagram of the age up-regulated genes (left) and age down-regulated genes identified between young and the old mice during sleep (over the four time points) in spontaneous sleep (SS) group or at ZT0 (FDR<1%). Figure S6. Heatmap of the log2 fold change between SD and SS for the 21 genes that are significant (FDR<5%) in the interaction tests (age*behavior state). Fold changes at each time point are shown for young (left) and old (center) mice, as well as the averaged fold changes (aveFC) for young and old mice across the four time points (right).The 12 sleep genes (top) and 9 wake genes (bottom) are further categorized into three classes, Y.specific, O.specific, or Common, based on if or not the genes are identified as sleep or wake genes in young animals, old animals, or both ages. Figure S8. Expression profiles of 12 core clock genes during spontaneous sleep (SS) and SD in the young (solid line) and old animals (dotted lines). Five genes significantly (FDR<0.01) altered in both young and old mice (Arntl, Per1, Per2, Dbp, and Tef) are marked by ‡. Five additional genes only significantly (FDR<0.01) altered in young animals (Clock, Cry1, Cry2, Hlf, and Nr1d2) are marked by ‡. By comparing young and old mice during SS or at ZT0, five genes showed overall significant age differences using moderated F-test (Cry1: FDR<0.01; Arntl, Per1, Per2, and Per3: FDR between 0.01-0.05), and asterisk (*) at individual time point indicate a nominally significant difference (P.value<0.05) between young and old. The multiple comparison corrected FDR values were provided in Table S3 for the moderated F-test of the SD vs SS contrasts over the four time-points for young and old mice, as well as the moderated F-tests of the O vs Y contrasts over the four time-points during SS and at ZT0. Figure S9. Expression profiles of genes related to synaptic homeostasis control during spontaneous sleep (SS) and SD in the young (solid lines) and in old animals (dotted lines). ‡ marks genes significantly altered between SD and SS in the young mice, and  marks genes significantly altered between SD and SS in the old mice. Although many transcription factors that are indicators of neuronal activities (Fos, Fosb, Egr1, Egr2, Arc, Bdnf, Ntrk2) are similarly regulated between sleep and wake in both young and old animals, many other genes playing key roles in synaptic assembly and neurotransmitter release were only regulated in the young but not in the old. Asterisk (*) at individual time point indicate a nominally significant difference (P.value<0.05) between young and old during SS or at ZT0. The multiple comparison corrected FDR values were provided in Supplement Table 4 for the moderated F-test of the SD vs SS contrasts over the four time-points for young and old mice, as well as the moderated F-tests of the O vs Y contrasts over the four time-points during SS and at ZT0. Mean±se fold change between SD and SS or between SS and ZT0 of the involved genes is shown for each enriched function. Fold changes were compared between young and old animals using paired Wilcoxon Signed Rank, and P-values < 0.05 are shown in bold italic. Effect sizes were calculated as Cohen's d indicating standardized difference between two means (of young and old mice). Effect sizes >=0.5 (medium effect) is shown in bold. *genes from three Signal/transmembrane functional clusters enriched among common, young-specific, and old-specific sleep genes were combined and compared between the age groups.