Reorganization of pancreas circadian transcriptome with aging

The evolutionarily conserved circadian system allows organisms to synchronize internal processes with 24-h cycling environmental timing cues, ensuring optimal adaptation. Like other organs, the pancreas function is under circadian control. Recent evidence suggests that aging by itself is associated with altered circadian homeostasis in different tissues which could affect the organ’s resiliency to aging-related pathologies. Pancreas pathologies of either endocrine or exocrine components are age-related. Whether pancreas circadian transcriptome output is affected by age is still unknown. To address this, here we profiled the impact of age on the pancreatic transcriptome over a full circadian cycle and elucidated a circadian transcriptome reorganization of pancreas by aging. Our study highlights gain of rhythms in the extrinsic cellular pathways in the aged pancreas and extends a potential role to fibroblast-associated mechanisms.


INTRODUCTION
Optimal biological functions of an organism are highly dependent on its ability to adapt to environmental pressures.Because of the rotation of the Earth around its own axis and the sun, life evolved under multiple rhythmic environmental regimes, including both seasonal and daily changes in light exposure [1].As a result, nearly every species on Earth has developed a biological timekeeping system that can anticipate these changes and organize physiology and behavior in a way that is advantageous for the organism.Biological rhythms of approximately 24 h, called circadian rhythms, are highly conserved among mammalian species and provide a temporal order for behavioral and physiological processes [2].Almost all physiological processes in mammalian organs are under circadian regulation [3,4].The circadian rhythm apparatus consists of a central master clock located in the suprachiasmatic nuclei (SCN) within the hypothalamus, which synchronizes peripheral oscillators located in various tissues such as the liver, pancreas, adipose tissue, and skeletal muscle.The SCN is entrained by light [5,6], which allows synchronization with the external environment i.e., the 24-hour light-dark cycle governed by the Earth's rotation.This directs the central and peripheral clocks to AGING adapt to changes in light, optimizing physiological processes to these daily cycles.Through several regulatory mechanisms (e.g., endocrine, neurological, thermal), the SCN coordinates responses with the peripheral clocks, which have their own phases, to ensure synchronized daily rhythms are maintained [7].In addition, peripheral rhythms can also be modulated, for example, by nutrient sensing (i.e., from food intake), hormonal cues and temperature8.Although the SCN acts as the master pacemaker in the human body, it is evident that circadian oscillations are observable in almost every cell of the body and these rhythms may persist in isolation from the SCN [8].
Studies have shown that pancreas function is strongly under circadian control and disruption of circadian rhythms exacerbates metabolic diseases that include Type II diabetes mellitus, obesity and metabolic syndrome in both animal models and humans [9][10][11][12].In addition, circadian functions are known to decline over the lifespan in several organs [13].For example, agerelated changes in the amplitude of sleep/wake activity, body temperature, metabolic functions, and hormone release is well documented [14].Aging is also connected with the inability to adjust to a new light/dark cycle, and increased mortality following repeated jet lag [15,16].A few recent age-related studies showed that the core clock in aged tissues remains largely intact, while the transcriptional output is significantly altered [17].
Given the strong connection between age and the resiliency of pancreas to pathologies, knowledge of differential pancreatic circadian transcriptome could lead to deeper understanding of its function and unearth novel targets for therapeutic interventions against agerelated pancreas pathologies such as metabolic disorders and cancers.Pancreas is composed of heterogenous cell types: most of the cells (90%) have exocrine function including production of digestive enzyme and the bicarbonate (by acinar and ductal cells, respectively), while endocrine cells comprise less than 2% of the pancreatic cells.The remainder of the pancreas mass is composed of the interstitial cells including various stromal and immune cells [18,19].Studies have shown an overall age-related decline in the function of pancreatic cells in conjunction with accumulation of cell-type specific transcriptional errors and oxidative stress [20].A systematic study of the age-related changes in the temporal gene expression profile of pancreas is largely lacking.
Here we carried out a 24-h circadian transcriptomic analysis of pancreas from male mice at young and old ages.We define a comprehensive circadian transcriptome landscape and identify biological pathways that are reflective of aging pancreas.
Additionally, analysis of pancreatic microenvironment revealed novel mechanistic insights into the fibroblast mediated regulation of rhythmicity in aged pancreas.We suggest that circadian transcriptome in aging pancreas re-organizes in response to age-specific signals from the cellular microenvironment, primarily modulated by fibroblasts.

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 2 > 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).

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 AGING 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.

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.

Fibroblast gene activity is rhythmic in aged pancreas
We hypothesized that in aged pancreas, fibroblasts in the microenvironment may be involved in reorganization 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 reorganization, 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 timestamped 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.

DISCUSSION
Despite evidence of strong circadian regulation of pancreatic functions as well as impact of aging on the circadian functions across several organs, the agedependent changes in temporal gene expression in pancreas remained unclear [9][10][11][12][14][15][16].To investigate this, we profiled transcriptome-wide gene expression of aging pancreas in mice and assessed the processes and mechanisms that gain or lose rhythmicity with age.
Our data showed that despite similar differential expression levels of steady state RNA (>98%), marked age specific differences exist in the circadian transcriptome of young and old pancreas (1,662 genes rhythmic), with extrinsic cellular processes gaining rhythm in the aged pancreas.These extrinsic pathways respond to maintain the architecture of the extracellular matrix (ECM), thus creating a three-dimensional network of macromolecules that provide structural and biomechanical support to tissue [24].Observed rhythmicity in ECM related processes specific to the aged pancreas indicates an age-specific re-organization of the transcriptome to ensure a temporal ECM integrity that could be pivotal towards pancreatic homeostasis and resilience in old age.
Interestingly, we found that the core circadian clock structure stayed largely conserved in aged pancreas, despite gains in rhythmicity of extrinsic cellular processes.This alludes to a local re-programming of the circadian transcriptome with aging in pancreas probably due to aging related stress response in the tissue [24].Further investigation into the-cause-and-effect connections between the changes in the ECM and rhythmicity of extrinsic cellular pathways could shed important insights into the mechanism of aging and could help identify novel targets to maintain healthy pancreas activity.
As a potential candidate that is involved in the ECM reprogramming, we looked at the cellular microenvironment encompassing immune and stromal cells (despite functional relevance to the endocrine function of pancreas, islet cells could not be investigated as they make up only 2% of the pancreatic mass 19 ).Fibroblasts were found to be the most abundant as well as the most rhythmic in aged pancreas, with direct gene-to-gene associations of rhythmicity to the extrinsic rhythmic pathways.This previously unknown temporal reorganization of pancreatic cellular microenvironment, centered on fibroblasts markers points to a fibroblastcentered circadian re-organization of the pancreas with aging.These findings are particularly interesting since fibroblasts are the major stromal population that expand in pancreatic adenocarcinoma, an aging-related pancreas pathology [28].Whether the circadian reordering of pancreatic microenvironment by aging could predispose pancreas to cancer by age, or whether similar circadian gain of the ECM and fibroblast signaling occurs in pancreatic cancer remains to be studied.

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 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 RNAcontaining 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.AGING 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).

Figure 1 .
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.

Figure 2 .
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.

Figure 3 .
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.