A systems genetics resource and analysis of sleep regulation in the mouse

Abstract

Sleep is essential for optimal brain functioning and health, but the biological substrates through which sleep delivers these beneficial effects remain largely unknown. We used a systems genetics approach in the BXD genetic reference population (GRP) of mice and assembled a comprehensive experimental knowledge base comprising a deep "sleep-wake" phenome, central and peripheral transcriptomes, and plasma metabolome data, collected under undisturbed baseline conditions and after sleep deprivation (SD). We present analytical tools to interactively interrogate the database, visualize the molecular networks altered by sleep loss, and prioritize candidate genes. We found that a one-time, short disruption of sleep already extensively reshaped the systems genetics landscape by altering 60%-78% of the transcriptomes and the metabolome, with numerous genetic loci affecting the magnitude and direction of change. Systems genetics integrative analyses drawing on all levels of organization imply α-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid (AMPA) receptor trafficking and fatty acid turnover as substrates of the negative effects of insufficient sleep. Our analyses demonstrate that genetic heterogeneity and the effects of insufficient sleep itself on the transcriptome and metabolome are far more widespread than previously reported.

Fig 1. Study design.

Thirty-three BXD lines plus the 2 parental strains and their reciprocal F1 progeny were phenotyped. Mice were submitted to either one of 2 experiments. In Experiment 1 (left), EEG/EMG signals and LMA were recorded under standard 12:12 h light–dark conditions (white and black bars under top-left panel) for 2 baseline days (B1, B2), a 6 h SD (red bar) from ZT0–6 (ZT0 = light onset), followed by 2 recovery days (R1, R2). The deep sleep-wake phenome consists of 341 sleep-wake state-, LMA-, and EEG-related phenotypes quantified in each mouse, among which time spent in NREM sleep (gray area spans mean maximum and minimum NREM sleep time among BXD lines, respectively, for consecutive 90 min intervals). Mice in Experiment 2 (right) were used to collect cortex, liver, and blood samples at ZT6. Half of the mice were challenged with an SD as in Experiment 1, the other half were left undisturbed and served as controls (labeled Ctr). Cortex and liver samples were used to quantify gene expression by RNA-seq, blood samples for a targeted analysis of 124 metabolites by LC/MS, or with FIA/MS. For phQTLs, mQTLs, and eQTLs, a high-density genotype dataset (Genome; approximately 11,000 SNPs) was created, merging identified RNA-seq variants with a publicly available database (www.genenetwork.org). The entirety of the multilevel dataset was integrated in a systems genetics analysis to chart molecular pathways underlying the many facets of sleep and the EEG, using newly developed computational tools to interactively visualize the results and pathways, and to prioritize candidate genes. EEG/EMG, electroencephalography/electromyogram; eQTL, expression quantitative trait locus; FIA/MS, flow injection analysis/mass spectrometry; LC/MS, liquid chromatography/mass spectrometry; LMA, locomotor activity; mQTL, metabolic quantitative trait locus; NREM, non-REM; phQTL, phenotypic quantitative trait locus; RNA-seq, RNA sequencing; SD, sleep deprivation; ZT, zeitgeber time.

Fig 2. Genetic diversity in the BXD panel greatly impacts behavioral, metabolic, and molecular traits.

The phenome was divided into 3 phenotypic categories: (i) LMA, (ii) EEG features (labeled EEG), and (iii) sleep-wake state characteristics (labeled State), which were subdivided further (see Materials and methods). The 5 classes of metabolites and the gene expression represent intermediate molecular phenotypic categories. (A) Heritability for EEG/behavioral and metabolite phenotypes. Dots represent single phenotypes within each category and subcategory indicated along the x-axis. Red dots represent phenotypes recorded in baseline (labeled bsl; B1 and B2), blue in recovery (labeled rec; R1 and R2), purple during SD, and green dots refer to the recovery-to-baseline contrasts. Values represent narrow-sense heritability. (B) Overview of significant and highly suggestive (FDR < 0.1) QTLs obtained for all 341 EEG/behavioral phenotypes (phQTLs: LMA in red, EEG in blue, and sleep-wake state in green) and 124 blood metabolite levels in baseline and recovery (mQTLs; purple). Note that overlap of neighboring QTLs renders color shading darker. (C) Venn diagram of genes under significant cis-eQTL effect in liver and cortex for the two experimental conditions (SD and controls [labeled Ctr]). EEG, electroencephalography; eQTL, expression quantitative trait locus; FDR, false discovery rate; LMA, locomotor activity; mQTL, metabolic quantitative trait locus; phQTL, phenotypic quantitative trait locus; QTL, quantitative trait locus; SD, sleep deprivation.

Fig 3. How to visualize multidimensional networks and prioritize candidate genes?

(A) Classical network visualization methods strongly depend on the layout algorithm used for positioning nodes, making structure interpretation and reproducibility difficult. (B) Hiveplot network visualization and structure strategy. See text for details. (C) The classical network visualization for the 3 phenotypes (blue nodes 1–3) in panel A can be represented with our method with 1 hiveplot per phenotype. Phenotype 1 showed more cortex–liver correlations than the 2 other phenotypes through 1 metabolite, connecting up- and down-regulated genes in cortex after SD and down-regulated genes in liver. Phenotype 2 shows genomic regions with strong allelic effect over multiple genes in liver and cortex through a high number of trans-eQTLs. Phenotype 3 was mostly connected to cortically expressed genes correlating strongly with up-regulated metabolites; most cis/trans-eQTLs affected only cortical genes. The number of significant (labeled sf) and suggestive (labeled sg) phQTLs detected for each phenotype are indicated on bottom left. The 3 phenotypes were related to active wake behaviors during recovery (Phenotype 1 and 2: LMA per hour awake and time in TDW, respectively, both during ZT12–24; Phenotype 3: Gain in time spent awake during ZT24–6). (D) Gene prioritization strategy to identify candidate genes associated with phenotype/metabolite variation, illustrated for 6 genes. Five types of analyses were integrated into a single score for each gene to reflect its strength as candidate gene, namely from left to right (i) and (ii) QTL mapping for gene expression (eQTLs) and ph- or mQTLs, respectively, (iii) DE after SD, (iv) gene expression/phenotype correlations, and (v) analysis of protein-damaging genetic variations relating genes to an allelic effect. See text for further details. (E) To illustrate and validate our scoring, strategy, genes in liver were prioritized for levels of α-AAA after SD. Dhtkd1 was identified as top-ranked candidate gene. Results from QTL mapping (red line) and prioritization analysis (green line); red and black horizontal lines indicate significant thresholds for the QTL and prioritization, respectively. α-AAA, alpha-aminoadipic acid; DE, differential expression; eQTL, expression quantitative trait locus; FDR, false discovery rate; LMA, locomotor activity; LOD, logarithm of odds ratio; mQTL, metabolic quantitative trait locus; QTL, quantitative trait locus; phQTL, phenotypic quantitative trait locus; SD, sleep deprivation; TDW, theta-dominated waking; ZT, zeitgeber time.

Fig 4. Profound effects of SD on transcriptome, metabolome, and phenome.

EEG/behavioral phenotypes, metabolites, and transcripts are organized into 3 “columns” (from left to right). Top 3 panels show the SD response (recovery/baseline fold change). Bottom 3 panels depict examples of allelic effects on the SD responses, with color-coding indicating the presence of a C57BL/6J or DBA/2J haplotype under the mapped QTL peaks (B6: gray for BXD and black for parental; D2: light brown for BXD and dark brown for parental). White bars mark the F1s and hatched bars strain in which haplotype could not be unambiguously determined. (A) Phenotypic changes after SD. The top significantly changed phenotype was the increase in NREM sleep EEG delta power (1–4 Hz) after SD (far-left blue data point). The most up-regulated phenotype was time spent in REM sleep during the first 6 h of darkness (ZT12–18) after SD (highest green data point). (B) Metabolite changes after SD. Most amino acids (blue) were down-regulated and most sphingolipids (brown) up-regulated after SD. The acylcarnitines C18:1 and C18:2 (highest red dots) increased the most. Vertical red line: significant threshold (FDR-adjusted p-value = 0.05). (C) DE analysis (SD/Ctr) for cortex (left) and liver (right). Genes were sorted according to their ranked p-value along the x-axis. Significantly affected transcripts in red (FDR-adjusted p-value < 0.05), nonsignificant results in black. Blue dots indicate 78 genes considered core molecular components of the sleep homeostatic response in the cortex [34]. Note that no low fold change threshold was applied. (D-F) Examples of genetically driven EEG/behavioral, metabolic, and transcriptional responses to SD, respectively. See text for details. Arc, activity-regulated cytoskeletal-associated protein; Ctr, control; DE, differential gene expression; EEG, electroencephalography; Egr2, early growth response 2; Fam107a, family with sequence similarity 107, A; FDR, false discovery rate; LMA, locomotor activity; Mlycd, malonyl-CoA decarboxylase; NREM, non-REM; Plin4, Perilipin 4; Pla2g4e, phospholipase A2, group IVE; QTL, quantitative trait locus; SD, sleep deprivation; Ttll8, tubulin tyrosine ligase-like family 8; ZT, zeitgeber time.

Fig 5. EEG delta power in NREM sleep after SD is associated with Kif16b and Wrn.

(A) NREM sleep EEG spectra in the first 3 h after SD (ZT6–9) for the 2 BXD lines that displayed the lowest and highest EEG activity in the fast delta frequency band (2.5–4.25 Hz, δ2; top, see panel E) and for the 2 BXD lines that displayed the smallest and largest increase (or gain) in EEG power in the slow delta band (1.0–2.25 Hz, δ1; bottom, see panel E). Spectra were “1/f-corrected” (and therefore not directly comparable to the values in panel E) for better visualization of activity in higher frequency bands (theta [5–9 Hz, θ], sigma [11–16 Hz, σ], beta [18–30 Hz, β], and slow [32–55 Hz, γ1] and fast gamma [55–80 Hz, γ2]). Subsequent analyses were performed without this correction. (B) QTL mapping and prioritization for δ2 power identified a significant association on chromosome 2 and Kif16b in cortex as top-ranked gene (top). For the δ1 increase after SD, we obtained a suggestive QTL on chromosome 8 and a significant prioritization score for the DNA-helicase Wrn. (C) Hiveplot visualization of network connections for the δ1 and δ2 power after SD (top-left panels) and the SD-induced increase in δ1 and δ2 power over baseline (bottom-left panels). Note the marked differences in the networks and QTLs regulating the expression of these 2 delta bands. Right hiveplots highlight Kif16b in the δ2 power–associated network (top), and Wrn in the network associated with the δ1 increase (bottom). Only Kif16b expression in the cortex was linked to the chromosome 2 cis-eQTL and was not associated with any metabolite. Wrn expression was significantly linked to the chromosome 8 cis-eQTL and to the long phosphatidylcholine, PC-ae-C38:5. (D) Kif16b is highly significantly down-regulated in cortex (left), while it remains unchanged in liver after SD (p = 0.15; not shown). Also, Wrn expression was strongly down-regulated by SD in cortex (right) and only marginally so, albeit significantly, in liver (p = 0.02; not shown). (E) Strain distribution patterns. BXD lines carrying a B6-allele on the chromosome 2–associated region showed higher δ2 power after SD (left) and a significantly higher Kif16b expression (p = 1.3e−15; second to left) than D2-allele carriers. D2-allele carriers of the chromosome 8–associated region showed a larger δ1 increase after SD (second to right) as well as a significantly larger decrease in Wrn expression after SD (right) than B6-allele carriers. For color-coding of genotypes, see Fig 4. CPM, counts per million; Ctr, control; EEG, electroencephalography; eQTL, expression quantitative trait locus; FDR, false discovery rate; Kif16b, Kinesin family member 16B; NREM, non-REM; PC-ae, phosphatidylcholine acyl-alkyl; QTL, quantitative trait locus; SD, sleep deprivation; Wrn, Werner syndrome RecQ like helicase; ZT, zeitgeber time

Fig 6. Changes in the frequency of theta oscillation during REM sleep after SD are associated with Cyp4a32.

(A) Spectral profiles of the REM sleep EEG for 2 strains displaying an opposite shift in the frequency of theta oscillations after SD relative to baseline. This shift was quantified by the decrease and increase in TPF for BXD61 and BXD101, respectively (see panel F). (B) Hiveplot for the SD-induced shift in TPF. (C) One significant QTL for the TPF shift was detected on chromosome 4 and 1 suggestive QTL on chromosome 8. Prioritization yielded Cyp4a32 as the top-ranked significant gene, based on the significant cis-eQTL modifying its expression in liver and a predicted damaging variation (V314E). (D) Effects of SD and genotype on liver Cyp4a32 expression. Carrying a B6-allele at the Cyp4a32 cis-eQTL–associated marker greatly decreased its expression. (E) Hiveplot for the SD-induced shift in TPF, highlighting Cyp4a32’s links to the amino acid Valine and the chromosome 4 eQTL marker. (F) Strain distribution patterns for TPF differences and liver Cyp4a32 expression after SD. B6-allele carriers at the chromosome 4–associated region had lower Cyp4a32 liver expression and a decrease in TPF after SD, while D2-carriers increase TPF and have higher Cyp4a32 expression. CPM, counts per million; Ctr, control; Cyp4a32, Cytochrome P450, family 4, subfamily a, polypeptide 32; DE, differential expression; EEG, electroencephalography; eQTL, expression quantitative trait locus; lod, logarithm of odds ratio; QTL, quantitative trait locus; SD, sleep deprivation; TPF, theta-peak frequency

Fig 7. NREM sleep gain in the first 6 h of the dark period after SD is associated with Acot11.

(A) Time course of hourly values of time spent in NREM sleep in baseline, SD (red area), and recovery for the 2 BXD lines showing the largest (BXD70; green) and lowest (BXD83; blue) NREM sleep gain during ZT12–18 (left). NREM sleep gain during 4 consecutive 6 h intervals during recovery compared to corresponding baseline intervals shows that in the recovery dark period (gray area), BXD83 mice did not accumulate extra NREM sleep, while BXD70 mice gained 88 min (middle). Strain distribution of ZT12–18 NREM sleep gain (right). B6-allele carriers compensated less for NREM sleep lost during SD than D2-allele carriers. For color-coding, see Fig 4. (B) Hiveplots for NREM sleep gain in 4 six-hour recovery intervals after the end of SD at ZT6. Compared to the other 3 intervals, NREM sleep gain was strongly associated with a number of metabolites during the second 6 h interval, i.e., ZT12–18. (C) NREM sleep gain during ZT12–18 mapped to a significant QTL on chromosome 4, explaining 45% of the total phenotypic variance (top left). PC-ae-C38:2 mapped suggestively to the same region (top right). Prioritization of liver transcripts for both phenotypes yielded Acot11 as top-ranked, significant gene (bottom). (D) Hiveplot for the ZT12–18 NREM sleep gain, highlighting Acot11. Acot11 was positively correlated with several phosphatidylcholines and to Ovgp1 expression in the cortex. (E) Allelic effect of the chromosome 4–associated region on Acot11 expression and PC-ae-C38:2 levels in the BXDs. Acot11 expression in liver after SD was under a strong eQTL effect (p = 1.6e−13) with B6-allele carriers showing a higher Acot11 expression than D2-allele carriers. B6-allele carriers also showed higher PC-ae-C38:2 levels after SD. (F) Both Acot11 and PC-ae-C38:2 levels changed after SD. Acot11 in liver and PC-ae-C38:2 in blood were significantly down-regulated. In the cortex, Acot11 was, however, significantly up-regulated, and the chromosome 4–associated region did not modulate cortical Acot11 expression. (G) Mice carrying 1 or 2 KO alleles for Acot11 displayed less extra NREM sleep during recovery. In contrast to the BXD panel, this difference was present in the second (ZT18–24, right) and not during the first (ZT12–18, left panel) 6 h of the recovery dark period. Acot11, acyl-CoA thioesterase 11; CPM, counts per million; Ctr, control; eQTL, expression quantitative trait locus; KO, knockout; NREM, non-REM; PC-ae, phosphatidylcholine acyl-alkyl; QTL, quantitative trait locus; SD, sleep deprivation; ZT, zeitgeber time