Reorganization of pancreas circadian transcriptome with aging

Changes in steady state gene expression in pancreas with aging

To quantitate the effects of aging on the molecular landscape of pancreas, we profiled the gene expression of the mice pancreas using RNA-seq (depth: >40 Mi reads). Two age-groups of male C57/BL6J mice were entered to our protocol (see method). Briefly, mice at 4 months or 24 months of age were singly housed and entrained under 12Light:12Dark conditions for 1 month. After the end of light phase of the last day of entrainment the lights remained off and the pancreatic tissues were collected in 4 hours intervals (6 time points, n = 2 mice per time point) starting at ZT0 in constant dark, at 5 months or 25 months of age. The transcriptome profiles for each time point were highly similar between age groups (R² > 0.95; Figure 1A). However, visualization over 24-h revealed clustering of young and old group separately, indicative of temporal molecular profiles (Figure 1B). Aged compared to young pancreas showed an upregulation in 133 genes and downregulation of 251 genes (p < 0.05; fold change ≥ 1.5) (Figure 1C and Supplementary Figure 1). Upregulated genes were enriched (FDR; p < 0.05) in pathways related to membrane transport and metabolism, consistent with the changes reported in the process of aging as well as the extracellular matrix (ECM) organization [21], while downregulated genes were enriched in pathways connected with the immune regulation and signaling (Rho GTPase, NF-κB) (Figure 1D).

Figure 1. Differential and temporal gene expression signatures in aged pancreas. (A) Pairwise correlation of normalized gene expression (transcripts per million; TPM) for young and old RNA-Seq samples (n=12). Colors indicate the Pearson correlation. (B) Principal Component Analysis (PCA) of temporal gene expression from young and old mice samples (n = 12). (C) Top: Volcano plot showing differential gene expression profile with significantly (p < 0.05, fold change ≥ 1.5) upregulated (n = 133 genes) and downregulated (n = 251) genes in old mice when compared to young. (D) Enriched molecular pathways corresponding to up and downregulated genes. (E) Top: Venn diagram of rhythmic genes in young and old mice. Bottom: Heatmap of z-scored rhythmically expressed genes in young and old mice. Heatmaps are double plotted for time.

Re-organization of circadian cycling of gene expression in aged pancreas

To delineate the differential circadian transcriptome of pancreas by age, we defined cycling genes of transcripts with 24-h cosine oscillations (Combined p-value; p < 0.05, Meta2D) [22]. We identified only 67 overlapping rhythmic genes, indicating that pancreas circadian transcriptome is age specific (Figure 1E). Notable rhythmic genes in both subgroups were the core clock genes, Arntl, Per2, Tef. On the other hand, genes involved in ECM reorganization, Col1a1 and Cdh11 were rhythmic only in old pancreas, while genes involved in mRNA metabolism were rhythmic in younger pancreas (e.g., Cpsf3, Figure 2A). Pathway enrichment analysis of rhythmic genes showed an enrichment in cell intrinsic processes (RNA metabolism, DNA replication) in young pancreas, while in the old group, the cell extrinsic cellular processes (ECM re-organization, collagen formation) were rhythmic (Figure 2B). Temporal pathway enrichment showed distinct processes enriched at different time points in young and old pancreas (Figure 2C), indicative of unique age specific circadian profiles. We then asked if this distinct circadian transcriptome could be due to age-dependent central circadian cues or molecular clock dysfunction. To assess the molecular clock, we first analyzed the hypothalamus, the location for the central circadian clock. The clock correlation distance (CCD) [23] of the temporal relative expression of 12 core clock genes revealed highly similar core clock profiles in both age groups (Supplementary Figure 2). An identical pattern also existed in the pancreas where clock CCD between young and old mice were similar (Figure 2D). These observations indicate that the core circadian clock progression is overall conserved, while extrinsic cellular processes gain rhythmicity in aged pancreas.

Figure 2. Circadian transcriptome of the aged pancreas. (A) Examples of circadian profiles of genes that are rhythmic in both young and old mice (Arntl, Per2, Tef), only in old mice (Col1a1, Cdh11), and only in young mice (Cpsf3). (B) Enriched molecular pathways corresponding to rhythmic genes in young and old mice. (C) Peak time maps of all oscillating genes from young and old mice. Each dot represents the peak time of 10 significantly rhythmic genes. Enriched pathways corresponding to each time point are also shown. (D) Clock correlation distance (CCD) of 12 core clock genes from young and old mice pancreas.

Fibroblast mediated age-modified cellular microenvironment is a hallmark of differentially rhythmic, aged pancreas

Extrinsic pathways are known to be modulated by the cellular microenvironment of interacting soluble factors and signals in the ECM consisting of stromal and immune cells [24]. We reasoned that the observed rhythmicity and upregulation in extrinsic pathways in aged pancreas could be a result of an age-modified tissue microenvironment. We imputed the gene signature based fractional abundances of 64 immune and stromal cell types from the old and young pancreas (Supplementary Figure 3A). Fibroblasts was estimated to be the most significantly altered cell type by far with markedly higher abundance in the aged pancreas (Figure 3A, p < 10^-8). Fluorescence Activated Cell Sorting (FACS) of the remaining pancreas of the studied mice conclusively showed significantly more fibroblasts in old mice (Figure 3A, CD90/Thy1+; p < 10^-9 and Supplementary Figure 3B). Additionally, CD90/Thy1, a marker of proliferative responses of fibroblasts to cytokines and growth factors [25], was found to be selectively rhythmic in old mice in our temporal gene expression analysis (Figure 3C). This data suggests that fibroblasts are heavily involved in the aging of pancreas, which is consistent with the known role of fibroblasts in aging [26] and points towards a fibroblast associated reorganization of pancreatic microenvironment. Using datasets of 178 human pancreas samples from Genotype-Tissue Expression (GTEx) project, we found a significant increase in fibroblast fractions in aged human pancreas (Figure 3B), indicating that this pattern is conserved across the species.

Figure 3. Fibroblasts mediate differential rhythmicity in aged pancreas. (A) Bar plots for (left) fibroblast fraction obtained by imputation and (right) Fibroblast to epithelial ratio obtained by Fluorescence Activated Cell Sorting (FACS) analysis for young and old mice pancreas. (B) Bar plots for fibroblast fraction obtained by imputation for young and old human pancreas from 178 GTeX datasets. (C) Heatmap of significantly rhythmic genes indicative of fibroblast markers in young and old mice. Fibroblast marker Thy1/CD90 is shown in bold. Heatmap represents the p-value of rhythmicity. (D) Top: Comparison of temporal fibroblast fractions obtained by two independent methods: imputation and FACS. Bottom: Similarity between the estimated fibroblast fractions between imputations and FACS methods. (E) Pairwise correlation plot for fibroblast markers and fibroblast responsive genes (FRGs) in young (lower diagonal) and old (upper diagonal) mice. Colors represent the Pearson correlation coefficient. (F) Heatmap of significantly rhythmic FRGs in young and old, time-stamped human datasets. Heatmap represents the p-value of rhythmicity. (G) Heatmap of significantly rhythmic FRGs in young and old mice. Heatmap represents the p-value of rhythmicity.

Fibroblast gene activity is rhythmic in aged pancreas

We hypothesized that in aged pancreas, fibroblasts in the microenvironment may be involved in re-organization of the rhythmic processes observed in Figure 2B. We assessed the temporal expression of genes reflective of fibroblast gene activity [27]. Interestingly, signatures for fibroblasts significantly encompassed rhythmic genes in aged pancreas (p < 10^-5, Figure 3C), with majority of hub proteins rhythmic including the fibroblast proliferation and activity marker CD90/Thy1 [25] (Supplementary Figures 3B, 4). No enrichment was observed in young pancreas. While there was a strong correlation between temporal fibroblast fractional estimates and FACS-based activity markers, they showed rhythmicity only in the old mice with similar 24-h cosine oscillations (Figure 3D).

Responsive genes to fibroblasts gain rhythm in the aged pancreas

Next, we speculated that genes that are responsive to fibroblast activity should show a similar, distinct pattern in the old pancreas. Using a panel of genes responsive to fibroblast gene activity signals [27] (Fibroblast Responsive Genes; FRGs, see methods), we found majority of FRGs depicted rhythmicity in aged pancreas (Figure 3G and Supplementary Figure 4). The rhythmic FRG enriched pathways were related to ECM re-organization, collagen formation, and integrin signaling. In addition, pairwise comparison of temporal expression profiles of fibroblast gene signatures and FRGs revealed a striking correlation (Figure 3E) suggesting rhythmic gene activity of FRGs mirrored fibroblast gene activity in aged pancreas. To test if a similar pattern exists in human pancreas, we assessed the rhythmicity of FRGs in time-stamped human pancreas datasets (from Ref 23). A significant number of FRGs were rhythmic in old but not in young human pancreas (Figure 3F). Overall, our study identified previously unknown circadian transcriptome reorganization of pancreas by aging, characterized by a gain of rhythm in the extrinsic cellular pathways linked to fibroblast-associated mechanisms.

Limitations of the study

Our work comes with some caveats that should be acknowledged. We collected gene expression data to assess rhythmicity over 1 cycle of 24h, at every 4h. This precluded repeated observation of peak-trough rhythmicity and could underpower our ability to detect more rhythmic genes. However, all the time points were collected in duplicates and mice were kept under constant darkness throughout the collection day to minimize the effect of light-induced changes.

In addition, we looked at any changes in rhythmicity profile across age groups in the form of pattern as reflected by the enrichment of biological pathways, and not an isolated single gene. This ensures that minor fluctuations in expression counts at a given time point for a given gene would not artificially alter the inference of rhythmicity. Lastly, further investigation into the connection between different cell populations of the cellular microenvironment and rhythmicity with aging is needed to obtain a comprehensive picture of the detailed mechanism by which ECM responsive processes gains rhythmicity in aged pancreas.

Animal care and use

C57BL/6J mice, born and raised at the University of Wisconsin-Parkside animal facility. The parental generation was obtained from the Jackson laboratory and the described procedures approved by the institutional animal care and use committee of the University of Wisconsin-Parkside. Breeding and housing colonies were housed under standard 12 hours of light followed by 12 hours of darkness (12L:12D) in light and temperature-controlled chambers and had continuous access to regular chow (ad-libitum) and fresh water. Mice were entered to our protocol at 4 months or 24 months of age. Both groups were singly housed and entrained for 4 weeks in LD chambers. Mice were sacrificed at 5 months or 25 months of age. Animal tissues following euthanasia were collected in 4-hour intervals under dark to allow for circadian analysis. Collected samples were partly persevered in RNA later for RNA processing, while a portion was used for flow cytometry as described below.

Sample preparation of total RNA sequencing

RNA was isolated from the pancreas using TRIzol extraction procedure. Briefly, frozen tissues were lightly ground in a mortar and pestle constantly submerged in liquid nitrogen. Frozen tissue between 10-100 mg was placed into RNase-free tubes containing TRIzol and RNase-free stainless-steel beads (NextAdvance, USA). Tissues were homogenized using a bullet blender at 4° C. Chloroform was used to separate the phases, and the RNA-containing aqueous phase was subjected to a column purification kit using RNeasy Mini RNA extraction kit (Qiagen, USA). RNA was DNase treated, and purified RNA was checked for integrity using an Agilent Bioanalyzer, all samples had an RNA integrity number above 8.0. RNAseq libraries were prepared using Illumina mRNA Prep kit (Illumina, USA). Libraries were pooled to equal molarity and sequenced on an Illumina NovaSeq (2 × 100bp) to achieve a minimum of 40M reads per sample. FastQ files were downloaded to the Rush University computing cluster for further data processing.

Flow cytometry

Pancreatic tissue was homogenized, and cells were suspended in DMEM and 10% fetal bovine serum (FBS). Supernatant was centrifuged at 500g for 5min, and the pellet was resuspended in a tissue digestion buffer (collagenase IV, dispase, DNase I, and DMEM). After arresting tissue digestion, the solution underwent centrifugation at 500g for 5min and the pellet was resuspended in DMEM, DNase 1, and 10% FBS. Total live cell retrieval varied based on timepoints, however ranged from 1.2-6.5 x 106 cells. Dead cells were excluded using the Live/Dead Violet Dead cell stain kit (dilution: 1:750; Thermo Fisher Scientific, USA). Cells were washed with 2% FBS and phosphate-buffered saline (PBS) and centrifuged at 2000rpm for 5pm. The pellet was resuspended in 100ul 2% FBS and PBS and incubated in Fc block for 5min in that dark at 4C to block nonspecific binding of antibodies. Cells then were incubated for 30 minutes at 4C in the dark with antibodies against cell surface markers. PE-conjugated anti-CD133 (dilution: 1:50; clone: 315-2C11), APC/Fire 700 anti- CD90.2- (dilution: 1:100; clone 5E10), FITC conjugated anti CD45 (dilution: 1:100; clone QA17A26), and APC/Fire 750 conjugated anti E-Cadherin (dilution: 1:100; clone DECMA-1. All antibodies purchase from BioLegend (USA). Samples were washed in 1mL 2% FBS and PBS and centrifuged at 2000rpm for 5mins. Pellets were resuspended in 400-500ul of cold adv DMEM, 10% FBS, and DNase 1, and the single cells were run on a BD LSRFortessa flow cytometer (BD Biosciences, USA). Data were analyzed using FlowJo software (Tree Star, USA).

RNA-seq data pre-processing

Reads from the RNA-Sequencing data were aligned to the Mus musculus genome assembly GRCm38 (mm10) using the STAR software [29]. Duplicated aligned reads were marked and removed using the Picardtools software. The gene expression count data were extracted using the HTseq software. The raw count data were normalized, followed by log2 transformation. We filtered out genes with mean read counts < 6. All data preprocessing was performed using R software.

Differential gene expression analysis

Genes with low counts were removed, and “expressed genes” were defined as those with at least 6 counts in 6 samples. VST was used to normalize RNA-seq across the 12 samples. PCA plots were used as an unbiased approach to visualize the distribution of clusters by age. DESeq2ref [30] was performed to address age-related changes in overall levels of gene expression in young vs. aged, and differentially expressed genes were defined as p-value < 0.05 and fold-change cutoff ≥ 1.5.

ANOVA

One-way ANOVA was calculated in R using car v.3.0.10ref [31]. Mean square difference between and within groups were calculated. Obtained F values were compared with the critical value in the F table to obtain P values. Inter-group differences were significant (p < 0.05) when the F value exceeded the critical F value for the given degrees of freedom.

Pathway enrichment

Pathway enrichment analysis for differentially expressed and rhythmic genes was performed with Reactome [32] using a hypergeometric statistical test and Benjamini-Hochberg FDR correction [33]. Redundant pathway terms were merged to create a parent term.

Circadian expression analysis

Meta2D cycling algorithm was used to search for genes whose expression cycle in a circadian manner8. Genes were considered cycling if the combined p values were 0.05 or lower for all three algorithms incorporated in Meta2D (ARSER, JTK_CYCLE and Lomb-Scargle).

Fractional cell abundance

Imputation of fractional abundances of 64 immune and stromal cell types was carried out using immunedeconv package in R. Immunedeconv provides a unified access to 10 most used computational methods for estimating immune cell fractions from bulk RNA sequencing data (quantiseq, cibersort, xcell etc.) [34].

Gene panels and datasets

Fibroblast gene panels were created from ref 27. Fibroblast responsive genes (FRGs) were identified by employing the fibroblast gene panel as input and annotated pathways in Reactome as matching criterion. The analysis was carried out in Cytoscape v.3.4.0. Hub proteins in a network were identified using CytoHubba (MCC, Neighborhood connectivity) plugin in Cytoscape.

GTeX datasets for human pancreas were obtained from GTeX Analysis V9 (dbGaP accession phs000424.v9). Time stamped human pancreas datasets were obtained from Ref 10 as transcripts per million (TPM).