Non-human primate seasonal transcriptome atlas reveals seasonal changes in physiology and diseases

Abstract

The metabolic, immune, and endocrine systems show profound seasonal changes in animals, including humans. In addition, morbidity from cardiovascular and psychiatric diseases is more severe and mortality rate is higher in winter. However, their molecular mechanisms remain unknown. Here we report the seasonal transcriptome of 80 tissues collected over 1 year from male and female rhesus macaques kept in a semi-natural outdoor environment. We find seasonal changes in plasma metabolites and hormones. Transcriptome analysis identifies sex differences in seasonally oscillating genes (SOGs) in all tissues studied, and we generate the web database 'Non-Human Primate Seasonal Transcriptome Atlas (NHPSTA).' Transcriptional regulatory network analysis, siRNA knockdown, and mutant mouse analyses reveal regulation of SOGs by GA-binding protein (GABP). We also demonstrate seasonal oscillations in the expression of disease risk factor genes and drug interacting genes. NHPSTA provides a molecular resource for seasonally regulated physiology and targets for therapeutic interventions for seasonally regulated diseases.

Fig. 1. Rhesus macaques show clear seasonal changes in physiology.

a Seasonal changes in testicular length in males and red coloration score for hindquarter skin in females. Blue and orange backgrounds represent winter and summer, respectively. b Heatmap showing seasonal changes in plasma metabolites measured by CE-MS. M: male, F: female. c Seasonal secretion profiles of plasma testosterone and estradiol. d Principal component analysis of transcriptomes of 924 samples across 80 neural and peripheral tissues. Colors indicate different tissue types. e Seasonal expression profiles of genes involved in the hypothalamus-pituitary-gonadal axis as determined by RNA-seq analysis. RS: the coefficient of determination.

Fig. 2. Identification of seasonally oscillating genes in male and female.

a Number of expressed genes in each tissue in both male and female b Number of SOGs in each tissue in both males and females. c Cumulative distribution of the peak phases of gene expression in different tissue types across 1 year.

Fig. 3. Comparison of seasonally oscillating genes between male and female macaques across 76 tissues.

Venn diagrams showing the numbers of SOGs from male (blue) and female (red) macaques.

Fig. 4. Seasonal transcriptions coordinate seasonal physiological function.

a GO enrichment analysis of common SOGs in both sexes. The top five enriched pathways in the following six tissues are shown: brown adipose tissue (BAT), skin on the back (SKNB), mesenteric lymph node (MEL), spleen (SPL), duodenum (DUO), and pars tuberalis of the pituitary gland (PT). b Heatmap showing the expression profiles of SOGs. Genes involved in representative pathways are highlighted.

Fig. 5. Regulation of SOGs by GABP.

a Distribution of SOGs ranked according to the number of tissues in which they show seasonal oscillation in male, female and both sexes. b GO enrichment analysis for the top 30 most common SOGs (one-sided Fisher’s exact test). c Distribution of potential transcription factors regulating SOGs ranked according to the number of tissues identified by transcriptional regulatory network analysis. d Effects of siRNA knockdown of GABPA and/or GABPB1 on GABPA and GABPB1 expression in macaque fibroblasts. (one-way ANOVA, Dunnett’s test, mean ± SEM, n = 4). e Effects of siRNA knockdown of GABPA and GABPB1 on the most common SOGs in macaque fibroblasts examined by qPCR (top) and RNA-seq (bottom) analyses. Numbers within the graph indicate the P value (unpaired two-tailed t-test, mean ± SEM, n = 4). f Comparison of 2,058 DEGs between GABPA and GABPB1 siRNA-treated and negative control siRNA-treated macaque fibroblasts (Benjamini-Hochberg correction adjusted P value (Padj) < 0.05 and fold change > 1.5). g Comparison of all 19,003 SOGs and 2058 GABP-regulated genes. h GO analysis of 1450 overlapping genes between SOGs and GABP-regulated genes (one-sided Fisher’s exact test).

Fig. 6. Seasonal tissue remodeling in mouse heart, pancreas, and kidney.

a Effects of 1-month exposure to short-day and cool (SC) or long-day and warm (LW) conditions on organ weight in CBA/N and C57BL/6 N mice. Numbers within the graph indicate the P value (unpaired two-tailed t-test, mean + SEM, n = 8). b (top) Effects of 1-week exposure to SC or LW conditions on Ki67-positive cell numbers in the heart (HEA), pancreas (PAN), and kidney (KID) of C57BL/N mice. Numbers within the graph indicate the P value (unpaired two-tailed t-test, mean + SEM, n = 4). (bottom) Effects of 1-week exposure to SC or LW conditions on cardiomyocyte diameter, pancreatic acinar cell size, and glomerular size. (t-test, mean + SEM, n = 4). c Representative images of Ki67 immunohistochemistry in the heart (HEA), pancreas (PAN), and kidney (KID) under SC and LW conditions. The number of Ki67 immuno-positive cells (red arrows) was higher under SC than LW conditions in the HEA and PAN, while glomerular size was larger under SC than LW conditions in the KID. d Representative images for the H&E staining of cardiomyocytes under SC and LW conditions. e Representative images of the immunohistochemistry of β-catenin, which localizes to the acinar cell membrane. Note that the nuclear density (i.e., number of nuclei in the field) is greater under LW than SC conditions, indicating a larger cytoplasmic area under SC conditions. L: Islet of Langerhans; G: Glomerulus. Scale bars: 50 μm. f Daily food and water intake in male C57BL/6 N mice kept under SC and LW conditions. Numbers within the graph indicate the P value (unpaired two-tailed t-test, mean + SEM, n = 6).

Fig. 7. In vivo role of GABP in mice.

a CRISPR/Cas9 targeting of the mouse Gabpb1 exon 9. Schematic representation of the mouse Gabpb1 gene illustrates crRNAs (blue characters) and the protospacer adjacent motif (PAM) sequences (red characters) tailored to exon 9. Sequencing of the deletion band shows ligation of the Cas9 target site. PCR across the genomic deletion region shows the smaller PCR product in mutants. Western blotting of Gabpβ1 with an antibody recognizing both long and short forms shows absence of long form of Gabpβ1 in the mutant homozygous (Mut/Mut) mice. β-actin was detected as loading control. b Effects of 1-week exposure to SC or LW conditions on Gabpb1 expression and organ weight in Gabpb1 exon 9 deletion mutant heterozygous mice and their wild-type littermates (two-way ANOVA, Šídák’s multiple comparisons test, mean + SEM, +/+: n = 3; +/Mut: n = 4). Numbers within the graph indicate the P value. c Immunohistochemistry of Ki67 in the heart (HEA), pancreas (PAN), and kidney (KID). In WT mice, the number of Ki67 immuno-positive cells (red arrows) was higher under SC than LW conditions in the HEA and PAN, while glomerular size was larger under SC than LW conditions in the KID. These differences were smaller in Gabpb1 heterozygous mutant mice. L: Islet of Langerhans; G: Glomerulus. d H&E staining of cardiomyocytes. Changes in cardiomyocyte diameter between SC and LW conditions were more obvious for WT than Gabpb1 heterozygous mice. e Immunohistochemistry of β-catenin, which localizes to the acinar cell membrane. Changes in the size of pancreatic acinar cells between SC and LW conditions were more obvious in WT than Gabpb1 heterozygous mice. Scale bars: 50 μm. f Effects of 1-week exposure to SC or LW conditions on Ki67-positive cell numbers, cardiomyocyte diameter, pancreatic acinar cell size, and glomerular size in Gabpb1 exon 9 deletion mutant heterozygous mice and their wild-type littermates (two-way ANOVA, Šídák’s multiple comparisons test, mean + SEM, +/+: n = 3; +/Mut: n = 4). Numbers within the graph indicate the P value.

Fig. 8. Seasonal changes in the expression of disease-related genes.

a DisGeNET enrichment analysis results for SOGs in the lung (LUN), aorta (AOR), bronchus (BRO), and colon (COL). b Heatmaps showing the expression profiles of SOGs. Genes involved in representative diseases are highlighted. c Seasonal expression profiles of ACE2 and TMPRSS2 in the BRO in males and females. d Heatmaps showing the expression profiles of SOGs. Risk factors for psychiatric diseases are highlighted.

Fig. 9. Seasonal changes in drug efficacy.

a Distribution of drugs ranked according to the number of targeted SOGs. All 19,003 SOGs were compared against the Drug Gene Interaction database. b Therapeutic categories of the top 50 drugs ranked by the number of targeted SOGs. c Effects of alcohol on spontaneous locomotor activity of mice kept under SC or LW conditions. Graphs depict the distance traveled in the open-field test at 1 and 3 h after EtOH administration. Numbers within the graph indicate the P value (two-way ANOVA, Šídák’s multiple comparisons test, mean + SEM, LW water treated group: n = 7, other groups: n = 8). d Effects of EtOH on the rotarod performance of mice kept under SC or LW conditions. Graphs depict the latency to fall from the rotarod as determined across five trials at 4 and 5 h after the administration. Numbers within the graph indicate the P value (two-way ANOVA, Šídák’s multiple comparisons test, mean + SEM, LW water treated group: n = 7, other groups: n = 8). e Effects of SC and LW conditions on EtOH metabolism. Numbers within the graph indicate the P value (two-way ANOVA, Šídák’s multiple comparisons test, mean + SEM, n = 8). f Distribution of SOGs ranked according to the number of targeting drugs.