The effects of graded levels of calorie restriction: VII. Topological rearrangement of hypothalamic aging networks

The effect of graded CR on the transcriptomic profile of aging-associated genes

We initially used an orthogonal signal correction partial least squares discriminant analysis (O-PLS-DA) to classify 12 and 24 hr AL fed groups (12AL, 24AL) and CR groups based on the expression levels of genes involved in aging, inflammation and oxidative stress (n = 408) obtained from a priori defined gene lists curated by Ingenuity Pathway Analysis (IPA) (‘inflammation of the nervous system’ and ‘oxidative stress response’) and GenAge database (‘aging-associated genes’ in Mus musculus). Individuals were clustered according to CR level, indicating that gene expression levels could predict CR level (Fig 1). Indeed model validation indicated that 41.8% of the variance in gene expression levels was explained by the different treatment groups (n = 6) (Xvar = 41.8 ± 18.5, p-value < 0.001). Of these 408 genes, expression analysis identified 16 genes differentially expressed at 24AL, 10 genes at 10CR, 14 genes at 20CR, 14 genes at 30CR and 59 genes at 40CR relative to 12AL (Table S1).

The effect of graded CR on network structure and connectivity

We estimated CMI based networks of these aging-associated genes and also standard gene co-expression networks based on correlations. CMI networks were characterized as having either a hierarchical or scale-free structure. Theory predicts that hierarchical networks should display a negative relationship between the degree to which genes cluster together (estimated as a cluster coefficient) and the number of neighbors of a given gene (node degree) (Fig 2A and 2B) [12]. The cluster coefficient did not decrease as node degree increased (Fig 2C), indicating that the networks were not hierarchical. Furthermore, node strength, a cumulative measure of node degree and the CMI of each gene, is predicted to follow a power distribution in scale-free networks (Fig 3A and 3B) [12]. However, in all networks node strength was not power law distributed (Fig 3C). In addition, visualizing the node strength on a log-log scale did not result in a linear relationship (Fig S1).

Another characterization of a network is the natural occurrence of gene clusters and the degree of connectivity of those clusters. This was estimated by modularity-based clusters of gene interactions and by assessing the strength of those interactions within these clusters. Modularity-based clusters of gene interactions changed and genes clustered differently at each CR level (Fig 4A) (Table S2). The modularity of the network, which represented the number of clusters, was also different at each CR level. The connectivity between genes in a cluster relative to the connections between clusters was measured by determining the modularity coefficient (Fig 4B). The highest modularity coefficient was observed at 10CR and 40CR, indicating that 10CR and 40CR had a higher connectivity between genes within clusters than between clusters (Fig 4C) (Table S2).

The effects of graded CR on gene centrality

We identified gene centrality, which are genes playing a central role in connectivity, by estimating gene eigenvector centrality in each CMI network. The variance in eigenvector centrality was reduced for 20CR and 30CR (Table S2) and we observed well-defined groups of genes with disproportionally large eigenvector centrality at 24AL, 12AL, 10CR and 40CR (Fig 5A) [35]. Central genes were identified based on the absolute maximum eigenvector centrality value (Table S2). At 24AL four genes were identified as central genes in the CMI network: excision repair cross-complementing rodent repair deficiency, complemen-tation group 4 (Ercc4), electron transferring flavoprotein, dehydrogenase (Etfdh), integrin alpha 9 (Itga9) and mastermind like 1 (Drosophila) (Maml1). One gene, ATP-binding cassette, sub-family C (CFTR/MRP), member 1 (Abcc1) was identified at 12AL, 10CR, and 30CR as a central gene. At 20CR three genes were identified as central genes: tumor necrosis factor (ligand) superfamily, member 10 (Tnfsf10), topoisomerase I binding, arginine/serine-rich (Topors), tripeptidyl peptidase II (Tpp2) and TNF receptor-associated factor 3 (Traf3). At 40CR the genes peroxisome proliferative activated receptor, gamma, coactivator 1 alpha (Ppargc1a) and palmitoyl-protein thioesterase 1 (Ppt1) were identified as central genes in the network. Ppargc1a and Ppt1 showed a clear graded increase in centrality, and Etfdh, Abbc1 and Traf3 a de-crease in centrality relative to CR level, but all genes did not show any simultaneous significant changes in gene expression levels relative to 12AL (Fig 5B). Furthermore at all levels of CR, genes with the highest eigenvector values were found in the same clusters (Fig 5C).

Comparison of networks based on CMI and correlation coefficients

We estimated Spearman's rank correlation coefficient between genes and constructed gene co-expression networks based on these correlations. The networks showed the same average number of edges between genes but their modularity and eigenvector centrality differed compared to CMI networks. Overall the gene co-expression network contained substantially more clusters compared to the CMI network with a higher modularity coefficient across CR. In contrast to the CMI based network, the modularity coefficient was the highest at 20CR (Table 1). The genes with the highest eigenvector centrality at 24AL, 20CR, 30CR and 40CR in the gene co-expression networks differed to CMI networks but were identical at 12AL and 10CR (Abcc1). Similar to the CMI networks, these genes were not significantly differentially expressed relative to 12AL. At 24AL angiotensinogen (serpin peptidase inhibitor, clade A, member 8) (Agt) was identified as central gene. In addition to Abcc1, four other genes were identified as central genes at 10CR: Agt, calcium/calmodulin-dependent protein kinase II gamma (Camk2g), CD200 antigen (Cd200), nuclear factor of activated T cells, cytoplasmic, calcineurin dependent 2 (Nfatc2) and transformation related protein 53 (Trp53). At 20 and 30CR the central gene was c-abl oncogene 1, non-receptor tyrosine kinase (Abl1). At 40CR the central gene was amyloid beta (A4) precursor protein (App).

Previous studies of gene expression networks in aging have mostly used gene co-expression networks to assess changes in network structure of individuals in relation to age and CR [16,22,23,36]. Here we used CMI because this measure provides a way to make inference about the amount of information shared between two genes when they are expressed. This is important for age-related studies for which we have a biological understanding that information flow at different cellular levels is disrupted by age. To our knowledge, CMI is used here for the first time to assess the changes in network topology associated with aging. We predicted that CR would result in a less fragmented regulatory network of aging associated genes compared to AL. We show that CR explained a significant proportion of the variability observed in the expression levels of aging-associated genes. Our CMI networks based on these genes suggest that (1) our networks were modular but did not have a hierarchical or scale-free topology; (2) the CMI network changed most at 40CR compared to 12AL with a higher connectivity between genes within clusters than between clusters; and (3) central genes for the regulatory network architecture showed a graded response to changing CR levels, with changes most pronounced at 40CR. This is consistent with previous studies at other biological levels for this experiment which showed a physiological state shift in mice exposed to 40CR [10,32–34,37]. Although the expression level of individual genes Ppargc1a, Ppt1, Etfdh, Abbc1 and Traf3 did not change [34], these genes were identified as ‘central genes’ indicating a network-wide role of these genes. We also showed that the networks exhibited a topological arrangement when exposed to graded CR.. It might be argued that when the networks are exposed to CR and close to a state shift, more non-monotonic relationships are formed which cannot be inferred using correlation methods.

Regulatory networks of aging-associated genes exposed to graded CR have a modular topology

Many biological networks have a scale-free topology displaying a high degree of clustering which is suggested to be a consequence of a hierarchical organization of that network [17]. The CMI networks of our mice were not scale-free and did not display a hierarchical topology. Instead we had a modular network characterized by highly interlinked modules (or clusters) which were connected to each other with relatively few links [38]. Both the number of clusters and the gene composition of clusters changed in response to CR. These changes in gene cluster membership would suggest changes in functional interactions between genes in response to CR [18]. In agreement with our results, a study with Caenorhabditis elegans found a structural reorganization of the transcriptomic network when glucose metabolism was impaired, and specifically in aging-relevant signaling pathways including the mammalian target of rapamycin mTOR pathway [39]. We also found that as CR levels increased the number of connections between genes was higher within gene clusters compared to those between clusters. The changes in connectivity observed here indicate that at 10CR and 40CR the networks have a higher connectivity in these functional gene groups compared to the two AL groups. The differences observed in network behavior at 20CR and 30CR could be a representation of a phase transition resulting in a topological state shift occurring at 40CR [40]. This would be in agreement with phenomic state shifts observed for these mice [33,37,86]. Prolonged “network stress” was found to induce a topological phase transition in cellular networks with an entirely different switch in cellular function (reviewed in [41]). Exposing mice to CR imposes a “network stress” and structural reorganization is necessary as an adaptive response. The adaptive response to CR is well established to increase longevity in many taxa and induce other beneficial anti-aging effects (reviewed in [42]). A “network stress” induced shift includes an initial decrease in average link density which might help the network to prevent propagation of damage [43]. During prolonged “network stress”, the network undergoes a segregation between less important genes and more important genes (central genes) which results in these central genes displaying more and stronger links to their neighbors. This would lead to more highly connected clusters of genes and less connections between genes outside these clusters which was observed at 40CR. As the level of CR increased the level of “network stress” increased and segregation occurred resulting in a reorganization and higher connectivity within clusters at 40CR. The network reorganization and increased connectivity at 40CR suggests the network would be more protected against the aging-observed decline in connectivity [22].

Changes in gene centrality irrespective of gene expression levels

Topologically central genes have a disproportionate influence on the overall integrity of gene regulatory networks [20]. Here, we identified central genes that played a central role in network structure and found a change in the connectivity of these genes in relation to the CR level. The central gene in the CMI networks at 12AL, 10CR and 30CR was Abcc1 or also known as MRP1. The importance of this gene in the network structure deceased in a graded manner with increasing CR level. Gene expression analysis of this gene showed no significant difference of this gene under CR compared to 12AL. This is in agreement with Laye et al. (2015) who demonstrated that the changes in correlation coefficient with age between metabolites in a co-expression network were independent of their mean levels [23]. Although first identified as a drug transporter, MRP1 is also involved in inflammation and oxidative stress [44]. MRP1 can also flux small peptides such as N-acetyl-leu-leu-norleucina (ANLL) which inhibits the 6S-proteasome activity which activates NF-ĸB [45]. Furthermore TNF-α has been associated with MRP1 [46]. We previously found that circulating TNF-α levels were reduced under CR [10] and that the transcription factor NF-ĸB in the hypothalamus was inhibited at all CR levels compared to AL, suggesting a reduced state of inflammation [34]. Furthermore micro-inflammation in the hypothalamus involving elevated NF-ĸB was recently suggested to mediate whole body-aging [5] and therefore a reduced state of inflammation under CR would have a protective role against whole-body aging.

At 20CR four genes were identified as central genes in network structure: Tnfsf10, Traf3, Tpp2 and Topors. The genes Tnfsf10 and Traf3 are associated with TNF-α and we established previously that Tnfsf10 correlated positively with circulating TNF-α levels [10,34]. Interestingly in human centenarians dementia is associated with high levels of circulating TNF-α levels and suggests a role of TNF-α in age-associated brain pathology [47]. Our results may indicate that reduced TNF-α signaling via Tnfsf10 and Traf3 plays a key role in network structure at 20CR. The other key gene, Tpp2, was shown to have a pro-longevity effect based on a mouse knockout study [48]. Tpp2 knockout mice over 1 year old exhibited elevated mortality associated with an aged appearance. Furthermore, alterations were observed in p53 expression and NF-ĸB activation but no direct molecular link was found with Tpp2 [48]. Tnfsf10 has also been found to mediate p53 dependent cell death [49]. Similar to Tpp2, Topors was also identified as having a pro-longevity effect [50] and also found to act as a tumor suppressor [51]. Topors knockout mice exhibited a significant reduction in mean and maximum lifespan and several mice had signs of premature aging [50].

The connectivity of the 40CR network differed substantially from the 12AL network. If our assertions are right, the 40CR network represents a network which receives the benefits of CR. Ppargc1a and Ppt1 were key topological central genes at 40CR which were involved in this overall difference in connectivity. They are important to protect the network from aging-associated loss in structure and fragmentation. A key change associated with aging was a loss in connectivity over the entire network and a decrease between genes within functional gene clusters [22]. Central genes are therefore expected to be particular foci for negative effects during the aging process [16]. Ppargc1a encodes PGC1α which is a major transcription factor contributing to gluconeogenesis, energy metabolism, lipid metabolism and functioning of the mitochondria [52]. PGC1α was shown to be upregulated in the brain when exposed to long term CR [53] but not in the hypothalamus when fasted for 60h [54]. Here, we did not find a significant regulation of Ppargc1a under CR compared to 12AL after 3 months of study. Brain aging was associated with a downregulation of the PGC1α mediated transcriptional pathway [55] and neuro-degenerative disorders such as Parkinson's disease are associated with lower levels of the target genes of PGC1α [56]. In Alzheimer's patients PGC1α expression is reduced [57]. Interestingly PGC1α overexpression in transgenic mice models for Alzheimer diseases exacerbates the neuropathological and behavioral deficits [58]. PGC1α null mice are lean and are resistant to diet-induced obesity. These mice exhibit profound hyperactivity associated with lesions in the brain region that controls movement [59]. When PGC1α was blocked locally within the central nervous system, fasted rats exhibited lower 24h food intake compared to those treated with a vehicle. Similar results were found for non-fasted high-fat induced obese rats but not for non-fasted low fat chow fed rats suggesting a greater role for neural PGC1α during fasting and obesity [60]. Fasted neural PGC1α null mice show diminished hypothalamic expression levels of the neuropeptides agouti-related neuropeptide (AgRP) and neuropeptide Y (NPY) which would suggest a role of PGC1α in regulating energy metabolism via regulation of hypothalamic neuropeptides [61]. The neuropeptides AgRP and NPY both increase hunger signaling and were upregulated in our transcriptomic dataset [34]. Although shown in diet induced obese mice, PGC1α is also associated with the control of reactive oxygen species (ROS) by hypothalamic pro-opiomelanocortin (POMC) neurons and hypothalamic ROS levels are positively correlated with circulating leptin levels [62]. Diano et al. (2011) postulated that the leptin resistance observed in these mice might be controlled by ROS manifested by reduced POMC and increased NPY/AgRP neuronal firing [62]. Interestingly during CR, leptin levels are reduced and correlated positively with reduced expression levels of POMC and negatively with increased expression levels of NPY/AgRP in our dataset [10,34] which suggests an improved leptin sensitivity and reduced levels of ROS. However, because of a lack of tissue we were unable to measure hypothalamic ROS levels in this study but we can conclude that our mice exhibited an increased ‘hunger profile’ [34]. Upon re-feeding, after a period of short-term (100 days) CR, the hyperphagic response suggests that hunger remained even after energy balance was re-established [63]. This elevated hunger profile might be a major factor contributing to the beneficial effect of CR [64]. A recent study has shown an abrogated effect on longevity when NPY^−/− mice were exposed to 30% CR compared to the AL fed mice, key role for NPY to link CR to increased lifespan [65].

We and others have previously postulated that hypothalamic nutrient sensing signaling might mediate the observed increase in lifespan under CR treatment [34,66–68]. Mice lacking PGC1α have abnormal diurnal rhythms of activity and metabolic rate and this was related to abnormalities in expression of clock genes and energy metabolism associate genes [69]. PGC1α also induces circadian clock resetting in restricted feeding and stimulates the clock gene aryl hydrocarbon receptor nuclear translocator-like (ARNTL) expression [69]. transcriptomic analysis and others indeed suggests that under CR mice are protected against aging-associated desynchronization of the circadian clock [34,36,70,71]. Although Ppargc1a gene expression levels were not significantly correlated with the extent of restriction, its downstream effects were. Here we were able to elucidate the importance of this gene in metabolic regulation via network analysis which would have remained unidentified with simple gene by gene expression based approaches. Interestingly mice with Ppt1 deficiency, the other key gene that changed in relation to the degree of CR, exhibited disruption of adaptive energy metabolism and downregulation of PGC1α [72]. Due to their key roles in network structure as central genes, Ppargc1a and Ppt1 play a major role in protecting the network from age-observed decline in connectivity [22,23].

Comparison of CMI and correlation coefficient based networks

Gene co-expression networks based on correlation coefficients have been widely used, for example in aging/CR research [16,23,36,73], to assess topological changes to networks in response to diseases [74–79], as characterization of a cellular state [80] and to understand plant cellular processes [81,82]. In contrast, CMI based networks have rarely been used despite recent theoretical advances in their application and interpretation and additional insight they can provide. Correlations assess ‘co-expression’ as the sum of the product of two random variables, while CMI specifically calculates the sum of the joint probabilities of two random variables from which inferences of causality can be made. Crucially, correlation assumes linear associations, or associations that change at the same rate with rank/order for non-parametric measures. CMI however, will detect associations whether they are monotonic or not. These measures are not antagonistic, they are complementary [83]. However, a major hurdle arises when we try to address the issue that in a many-body system the association between any given pair of objects can be influenced by their own associations with other objects. In our case this represents our attempt to detect whether the co-expression of two genes is real or caused by their mutual co-expression with a third gene. CMI estimates are more robust than partial correlation estimates because we have a better understanding of how to define joint probability distributions. Therefore, CMI can be more advantageous, particularly for large networks [84]. Finally correlation and CMI capture the shared information/covariance between two different measures of gene expression. For correlation, we try to understand the shared variance in gene expression. With CMI we try to estimate the joint information emanating from genes and thereby assuming causality. These two methods have previously been compared and an almost one-to-one correspondence between the two methods was found, the authors concluding that non-linear relationships did not play an important role in their gene networks [83,85]. The fundamental difference however lies in what can be inferred from these two methods. We indeed found similarities between our CMI and gene co-expression networks with respect to eigenvector centrality at 12AL and 10CR but we also found fundamental differences in modularity. In the CMI network the intermediate groups had fewer clusters but the amount of clusters at 40CR was more similar to 12AL while modularity differed. This was clearly shown by the modularity coefficient of the CMI based analysis while this was not the case for the correlation coefficient based networks. We can conclude from these results that gene co-expression based modularity coefficient estimates were relatively insensitive to CR whereas the CMI estimates were able to discriminate between 20CR and 30CR. It might be argued that when the networks are exposed to CR and close to a state shift, more non-monotonic relationships are formed which cannot be inferred using correlation methods. Hence CMI provides a measure for which we can directly test biologically relevant co-expression and given the particular focus we can measure changes in information flow and regulation with aging and CR.

In this study we found a structural reorganization of the gene regulatory network under graded levels of CR. Using the CMI network analysis approach we identified Ppargc1a, Ppargc1a, Ppt1, Etfdh, Traf3 and Abcc1 as key regulatory genes in network structure with changes in their centrality as the level of CR increased. Since expression levels of these genes did not change with CR, this study emphasizes the importance to move beyond gene by gene level analyses to better understand regulatory changes of the transcriptome when it is challenged (in our case by CR). Our results suggest that the network structure of aging-associated genes under graded CR was altered and might play an important role in preventing the loss of network structure observed with age.

Animals and experimental manipulations

All procedures were approved by the University of Aberdeen ethical approval committee and carried out under the Animals (Scientific Procedures) Act 1986 Home Office license (PPL 60/4366 held by JRS). Forty nine male C57BL/6 mice (Mus musculus), purchased from Charles River (Ormiston, UK) were individually housed and free access to water was provided. Mice were exposed to 12 hour dark/light cycle (lights on at 0630h) and body mass and food intake were recorded daily, immediately prior to nocturnal feeding. At 20 weeks of age (resembling early adulthood in humans), mice had their food intake monitored for 2 weeks of baseline and were then randomly allocated into 6 different treatment groups: 24h ad libitum intake (24AL) (n=8), 12AL intake which had access to ad libitum food but for only for the 12 hours of darkness per day, sometimes also called ‘time restricted feeding’ (n=8), 10 CR (n=8), 20CR (n=8), 30CR (n=8) and 40CR (n=9). Mice exposed to 40CR indicates 40% fewer calories were provided than the individual intake measured over a baseline period of 14 days prior to introducing CR. This is a CR protocol rather than a caloric dilution experiment [87].

Animals fed completely AL (i.e., having 24 hours access to food) may potentially over feed, become obese and CR associated changes compared to 24AL would therefore most likely reflect the anti-obesity effect of CR [42,88]. To address this issue, 12AL was set as a reference and graded levels of CR were introduced to investigate a potential graded response mirroring the graded lifespan response to CR manipulation [87,89]. Detailed information on overall study design, diet composition and detailed rationale are described elsewhere [32].

RNA isolation, RNA-sequencing, alignment and analytical procedure

After culling by a terminal CO₂ overdose, brains were removed, weighed and frozen in isopentane over dry ice and stored at −80°C for RNA isolation. The hypothalamus was carefully dissected at a later stage and RNA was isolated by homogenizing in Trizol (Sigma Aldrich, UK) according to manufacturer's instructions. Prior to RNA quantification, using the Agilent RNA 6000 Nano Kit, samples were denatured at 65°C.

Due to the very small size of the hypothalamus, some samples did not contain sufficient quantity of high quality RNA. In total, the RNA of 37 individual mice (12 h AL n=6, 24 h AL n=6, 10 % CR n=7, 20 % CR n=5, 30 % CR n=5, 40 % CR n=8) was successfully isolated and sent to Beijing Genomic Institute (BGI, Hong Kong) for RNA sequencing. Library preparation was done using a standard protocol of BGI and the library products were sequenced using an Illumina Hi-seq 2000, resulting in 50 bp single end reads. Standard primers and barcodes developed by BGI were used. Detailed information on library preparation for these samples and the bio-informatic pipeline has been described previously [34].

Network construction and biological interpretation

Based on the prior knowledge that aged brains have increased oxidative stress and inflammation [2,3], we curated a data set based on a prior defined gene lists from Ingenuity pathway analysis (IPA) (Ingenuity® Systems, http://www.ingenuity.com/products/ipa,version 2000-2015): ‘inflammation of the nervous system’ (n=196) and ‘oxidative stress response’ (n=155); and genes associated with aging derived from the online database GenAge (build 18, version 11/10/2015) (n=99) [48]. Normalized read counts (for method see [34]) of these a priori identified gene were initially analyzed by an O-PLS-DA. The Devium package remove the first package with PLS and OPLS R command functions were used (https://github.com/dgrapov/devium.git) to validate the model in the statistical environment R (version 3.1.2) [90].

The Bioconductor package minet [91] was used to infer large transcriptomic networks using CMI. The package returns a network where the nodes are genes and the edges are statistical dependencies between genes (CMI matrix). Four different entropy estimators and four different inference methods are available of which we used the Miller-Madow asymptotic bias corrected empirical estimator [92] and the ARACNE method [26]. CMI between genes was inferred in two steps: (1) MI is estimated between all pairs of genes and (2) the ARACNE inference algorithm applies a filter to remove estimated spurious links. ARACNE initially assigns to each pair of genes (i.e. nodes) a weight equal to MI and then removes any indirect interaction [26]. A network based on Spearman correlation coefficients between gene expression levels were also estimated with minet.

Gene centrality measures were estimated for each network based on eigenvector values. Eigenvector centrality captures the relative importance of each gene for the overall interaction network structure by estimating the contribution of the variance in interaction of each gene to the overall interaction heterogeneity in the entire network [21]. Gene eigenvector centrality was estimated by Eigen-decomposing each CMI network and using the elements of the eigenvector corresponding to the dominant eigenvalue from these decompositions. Gene strength per gene was calculated by the sum of the CMI between that gene and other genes in the network (row sums of the CMI matrix).

Weighted clustering coefficient estimates were used as described in [93]. Clusters of genes were identified based on defining a parsimonious division of a network [94,95]. This approach maximizes the number of connections (edges) within a cluster and minimizes the number of connections between clusters. The modularity coefficient quantifies the modularity technique and is a sum of associations for all genes belonging to the same clusters minus its expected value if genes were associated at random, based on the strength of each gene (method based on [93]).