Murine Gut Microbiome Meta-analysis Reveals Alterations in Carbohydrate Metabolism in Response to Aging

ABSTRACT Compositional and functional alterations to the gut microbiota during aging are hypothesized to potentially impact our health. Thus, determining aging-specific gut microbiome alterations is critical for developing microbiome-based strategies to improve health and promote longevity in the elderly. In this study, we performed a meta-analysis of publicly available 16S rRNA gene sequencing data from studies investigating the effect of aging on the gut microbiome in mice. Aging reproducibly increased gut microbial alpha diversity and shifted the microbial community structure in mice. We applied the bioinformatic tool PICRUSt2 to predict microbial metagenome function and established a random forest classifier to differentiate between microbial communities from young and old hosts and to identify aging-specific metabolic features. In independent validation data sets, this classifier achieved an area under the receiver operating characteristic curve (AUC) of 0.75 to 0.97 in differentiating microbiomes from young and old hosts. We found that 50% of the most important predicted aging-specific metabolic features were involved in carbohydrate metabolism. Furthermore, fecal short-chain fatty acid (SCFA) concentrations were significantly decreased in old mice, and the expression of the SCFA receptor Gpr41 in the colon was significantly correlated with the relative abundances of gut microbes and microbial carbohydrate metabolic pathways. In conclusion, this study identified aging-specific alterations in the composition and function of the gut microbiome and revealed a potential relationship between aging, microbial carbohydrate metabolism, fecal SCFA, and colonic Gpr41 expression. IMPORTANCE Aging-associated microbial alteration is hypothesized to play an important role in host health and longevity. However, investigations regarding specific gut microbes or microbial functional alterations associated with aging have had inconsistent results. We performed a meta-analysis across 5 independent studies to investigate the effect of aging on the gut microbiome in mice. Our analysis revealed that aging increased gut microbial alpha diversity and shifted the microbial community structure. To determine if we could reliably differentiate the gut microbiomes from young and old hosts, we established a random forest classifier based on predicted metagenome function and validated its performance against independent data sets. Alterations in microbial carbohydrate metabolism and decreased fecal short-chain fatty acid (SCFA) concentrations were key features of aging and correlated with host colonic expression of the SCFA receptor Gpr41. This study advances our understanding of the impact of aging on the gut microbiome and proposes a hypothesis that alterations in gut microbiota-derived SCFA-host GPR41 signaling are a feature of aging.

IMPORTANCE Aging-associated microbial alteration is hypothesized to play an important role in host health and longevity. However, investigations regarding specific gut microbes or microbial functional alterations associated with aging have had inconsistent results. We performed a meta-analysis across 5 independent studies to investigate the effect of aging on the gut microbiome in mice. Our analysis revealed that aging increased gut microbial alpha diversity and shifted the microbial community structure. To determine if we could reliably differentiate the gut microbiomes from young and old hosts, we established a random forest classifier based on predicted metagenome function and validated its performance against independent data sets. Alterations in microbial carbohydrate metabolism and decreased fecal short-chain fatty acid (SCFA) concentrations were key features of aging and correlated with host colonic expression of the SCFA receptor Gpr41. This study advances our understanding of the impact of aging on the gut microbiome and proposes a hypothesis that alterations in gut microbiota-derived SCFA-host GPR41 signaling are a feature of aging.
T he human gut harbors a complex ecosystem of microorganisms, collectively referred to as the gut microbiota. The initial colonization of the gut microbiota is generally thought to be established at birth, developing toward a stable adult-like microbial community structure by the age of 3 years (1). The gut microbiota develops with human hosts and has a variety of biological functions in human health and disease. Increasing evidence suggests that the gut microbiota plays important roles in the development and maturation of the immune system, food digestion and nutrient metabolism, gut endocrine regulation, neurological signaling, defense against pathogens, and elimination of toxins (2).
Aging is a time-dependent multifactorial process characterized by functional declines in a variety of biological and physiological systems and an increased risk for major human chronic diseases (3). During the aging process, remarkable aging-related gut microbial alterations have been observed in human studies across various populations (4,5). However, gut microbiome profiles across different elderly populations vary widely. A number of studies have found that the aged human gut microbiome is characterized by increased proportions of opportunistic pathogens (e.g., Clostridioides) and decreases in health-promoting bacteria (e.g., bifidobacteria and lactobacilli) (4)(5)(6)(7)(8). Moreover, several studies conducted on extreme-aged people (centenarians) show that the microbial profiles of extreme-aged people differ from those of old adults, potentially identifying a healthy aging microbial pattern (9)(10)(11). This healthy aging microbiome is characterized by the depletion of core microbial taxa (e.g., Bacteroides and Prevotella spp.) and becomes increasingly unique to individuals (9). The heterogeneity of findings in human aging studies indicates that multiple gut microbiome patterns of aging exist.
As we age, our lifestyle and diet are likely to change due to a reduction in our physiological systems. Aging is often accompanied by reduced food diversity, especially in fiber-containing foods (e.g., fruits and vegetables), as well as an increased risk of malnutrition. A healthy and diverse diet is positively related to a more diverse gut microbiota in old people (12). Dietary supplementation of prebiotics (e.g., inulin) leads to a higher abundance of beneficial bacterial groups (e.g., Bifidobacterium) in old people (13). On the other hand, malnutrition is associated with an increased abundance of Clostridiales taxa, which has been associated with increased frailty in aged individuals (12,14). Thus, diet may account for at least a portion of the variation in microbial composition among different elderly populations. Moreover, aged populations are more likely to experience hospitalization, infections, and exposure to antibiotics and other medications (15), which would also significantly impact the gut microbiome composition. The gut microbiota profile of old people in long-term-care facilities differs from those of community-dwelling subjects, with a decrease in microbial diversity and a loss of community-associated taxa (e.g., Coprococcus and Roseburia) (12). Thus, the interpretation of studies correlating age with changes in microbial communities is complex due to the various factors that would potentially impact and confound our understanding of aging-associated gut microbiome alterations. Systemically investigating and separating the effects of aging and exposures on gut microbiota composition and function are difficult in humans but may be accomplished in experimental animal models.
Research animals, such as mice (Mus musculus), are maintained under standardized conditions, with uniform diet, bedding, dark-light cycles, and other environmental exposures. Thus, experimental mice are a useful model system to study aging-specific changes in the host-microbiome interaction. Aging-related gut microbial alterations have been investigated in mice in several studies (16)(17)(18)(19)(20). A distinct gut microbial composition and a structure shift between young and old mice have been found, which partially resemble the changes observed in human studies. However, when we examined changes in specific taxa or functional alterations, inconsistent or conflicting results were common. For example, the Firmicutes-to-Bacteroidetes ratio (F/B ratio), an index often used to describe gut microbiome profiles, was increased in aged mice in a study by Kim et al. (21), while it was decreased in a study by Zhang et al. (22) and unchanged in a study by Langille et al. (23).
Inconsistent findings from different mouse studies make it difficult to come to a consensus regarding aging-specific microbial alterations, a prerequisite to identifying microbiome-based therapeutic targets for aging-related diseases. In this study, we performed a meta-analysis of publicly available 16S rRNA gene sequencing data of studies investigating the effect of aging on the gut microbiome in mice. By integrating agingrelated microbiome alterations across 5 independent studies, we hoped to identify age-specific changes in community composition or function common to all studies and then validate these in additional data sets. In brief, we reanalyzed raw sequencing data from all studies using the same pipeline. We then applied the bioinformatic tool PICRUSt2 (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States) to predict microbial metagenome function and established a random forest classifier to differentiate young and old age and identify aging-specific metabolic features. We then validated our random forest classifier against one 16S rRNA gene sequencing data set and one whole-genome metagenomic sequencing data set from young and old mice generated by our laboratory. Our analysis identified an increase in gut microbial diversity with age and revealed aging-associated alterations in microbial carbohydrate metabolism that significantly correlated with colonic Gpr41 expression, revealing how aging affects one mechanism of host-microbe interactions, that is, gut microbiome-produced short-chain fatty acids (SCFA) acting through their receptors.

RESULTS
Microbial alpha diversity increased with aging. The bioinformatic analysis pipeline used to process and analyze the raw sequencing data is shown in Fig. 1 and described in detail in Materials and Methods. Common metrics for alpha diversity, including Shannon, Simpson, Faith's phylogenetic diversity (PD), and Chao1 indices, FIG 1 Bioinformatic analysis pipeline. The raw 16S rRNA gene sequencing data were downloaded and imported into QIIME2 for data processing. DADA2 was used to denoise reads. OTU picking was performed by the Vsearch closed-reference method. The resulting feature table and feature data were then used for phylogenic tree construction and taxonomy annotation, followed by diversity analysis and taxon bar plot visualization. Functional prediction was performed by PICRUSt2. A random forest classifier was established based on predicted functional pathways to differentiate young-and old-age gut microbiomes and identify aging-specific metabolic features. The random forest classifier was then verified by two external data sets generated by us. The CB6F1 16S rRNA gene sequencing data were analyzed according to the same analysis pipeline. The whole-genome sequencing data were first filtered to remove low-quality and host genome contaminant reads, followed by SortMeRNA and OTU picking for taxonomy characterization and HUMAnN3 for functional analysis.
Aging Murine Gut Microbiome Meta-analysis mSystems were calculated. Due to the heterogeneity of studies, including varied sample processing, sequencing platforms, variable regions, study sizes, and sequencing depths across data sets (detailed in Table 1 and Fig. S2 in the supplemental material), the data were visualized by log 2 fold changes (log 2 FCs) of old to young with 95% confidence intervals (CIs) on a per-data-set basis (Fig. 2). With the exception of experiment 5 (Exp5), each data set had one or more alpha diversity indices that were significantly increased in the aged mice (Fig. 2). When considering all data sets together, all 4 alpha diversity indices were significantly increased in the aged mice ( Fig. 2) (24). Contradictory results regarding the direction of F/B ratio alterations in response to aging have been reported (21)(22)(23). Thus, we calculated the relative abundance of taxa at the phylum level and the F/B ratio. As expected,  Firmicutes and Bacteroidetes dominated the microbial community in both young and aged mice (Fig. 3A). However, no alteration of the F/B ratio was found in the aged mice compared to the young mice when considering all the data in aggregate (Fig. 3B)   Aging shifts gut microbiome structure. The beta diversity indices unweighted and weighted UniFrac distances were calculated and visualized with principal-coordinate analysis (PCoA) plots per data set and for all data sets in aggregation ( Fig. 4 and Fig. S3). Permutational multivariate analysis of variance (PERMANOVA) with 999 permutations was used to test for significant differences between young and old mice. Based on unweighted UniFrac distance matrices, all data sets demonstrated a significant shift in the gut microbiome between the young and old mice (P , 0.05) (Fig. 4A). When all the data sets were aggregated, a significant age-dependent alteration of the microbial community structure was observed based on unweighted UniFrac distance matrices (P = 0.001) (Fig. 4B). However, based on weighted UniFrac distance matrices, no significant difference was observed in 2 of 6 data sets (Fig. S3) or in the aggregate data (P = 0.259) (Fig. 4C). These results indicate that rare taxa rather than dominant taxa are primarily responsible for the microbial structure differences between young and old mice.
Altered gut microbial carbohydrate metabolism in response to aging. Random forest was used to differentiate young and old age based on operational taxonomic units (OTUs) and predicted metagenomic function by PICRUSt2 to identify aging-specific microbial signatures. An iterative random forest (iRF) algorithm was performed to identify features that were stable across 10 iterations, of which the iteration with the smallest out-of-bag (OOB) error rate was chosen as the final classifier for further analysis. We calculated the area under the receiver operating characteristic (ROC) curve (AUC) to estimate the performance of the random forest classifier. The AUC values of the model based on predicted metagenomic function for the test data set, the external CB6F1 16S data set, and the external B6 metagenomics data set were 0.75, 0.97, and 0.79, respectively ( Fig. 5A and B). Interestingly, the model based on OTUs did not perform as well as the one based on predicted metagenomic function, which achieved  S5). This result could be explained by microbial functional redundancy (25). That is, despite large variations in taxonomic composition across studies of mice from different vivariums, consistent changes in functional capacity in response to aging can be observed. Thus, metagenomic functional features selected by the random forest model might be more reliable than taxonomic features to stratify young and old mice and exhibit similar alterations in response to aging across studies. Next, aging-specific metagenomic features were selected based on the ranked mean decrease in Gini impurities (mean decrease Gini) (Fig. 5C), which measures the importance of a feature for classification in the random forest model. Interestingly, we found that 50% of the 20 most important features are involved in carbohydrate metabolism, including mannan degradation, gluconeogenesis I, and glucose and glucose-1-phosphate degradation. We also performed an enrichment analysis on predicted metagenomic function pathways across the studies (linear mixed-effect model adjusted for study). Despite the different methods performed, we found that 13 out of 25 significantly enriched pathways were involved in carbohydrate metabolism (Fig. S6) (BH-corrected P , 0.05). Taken together, alterations in the carbohydrate metabolism potential of the gut microbial metagenome are a key distinguishing feature of aging.
Decreased fecal SCFA in response to aging. Given the altered microbial carbohydrate metabolism metagenome with aging described above, we investigated whether fecal short-chain fatty acids (SCFA), the major end products of microbial carbohydrate fermentation, were changed with aging as well. We measured three primary SCFA concentrations in fecal samples from young and old mice used to generate the metagenomic data set. Fecal acetate (P , 0.001), propionate (P , 0.01), butyrate (P , 0.001), and total SCFA (P , 0.001) concentrations were significantly decreased with age (Fig. 6A). The decreased SCFA concentrations with aging in our cohort are consistent with previous studies in humans, mice, and rats (26)(27)(28)(29). Additionally, Spearman correlation analysis was performed between fecal SCFA and taxa at the genus level (Fig. 6B). We found several taxa considered to be SCFA producers or beneficial microbes (30)(31)(32)(33) that were positively correlated with one or more fecal SCFA, including Dorea (BH-corrected P , 0.05 for acetate, propionate, butyrate, and total SCFA), Oscillibacter (BH-corrected P , 0.05 for acetate, butyrate, and total SCFA), Mucispirillum (BH-corrected P , 0.05 for butyrate), Acetatifactor (BH-corrected P , 0.05 for butyrate and total SCFA), and Anaerotruncus (BH-corrected P , 0.05 for butyrate and total SCFA). Additionally, several taxa were significantly anticorrelated with one or more fecal SCFA, including Bifidobacterium, Alistipes, Parabacteroides, Prevotella, Dubosiella, Turicimonas, Bacteroides, Blautia, Erysipelatoclostridium, and Romboutsia (BHcorrected P , 0.05 for one or more metabolites). The negative correlation between Bifidobacterium and fecal SCFA is unexpected as Bifidobacterium is well known as a beneficial microbe that metabolizes dietary fibers to acetate and lactate via the bifid shunt (34). Bifidobacterium may alter metabolic pathways in response to environmental conditions (e.g., the availability of carbohydrates) (35). For example, Bifidobacterium could shift toward formate production at the expense of lactate when there is an increased need for ATP (35,36). Furthermore, a correlation does not necessarily indicate a direct relationship between a given microbe and SCFA generation. Further studies are needed to investigate the direct contribution of particular microbes to fecal SCFA via in vitro microbial community fermentation or investigations in vivo utilizing gnotobiotic mice.
Correlation of SCFA receptor Gpr41 gene expression with carbohydrate metabolism pathways. G-protein-coupled receptor 41 (GPR41) and GPR43 are the primary receptors that sense SCFA and mediate the biological effect of SCFA on host health (37). To assess whether the decreased fecal SCFA with aging would impact the expression level of SCFA receptors, we measured Gpr41 and Gpr43 gene expression in the colons of young and aged mice. Gpr41 expression was significantly decreased while Gpr43 expression was significantly increased in the aged mice compared to young mice (P , 0.05) (Fig. 7A). Gpr41 expression was significantly correlated with carbohydrate degradation and fermentation pathways, including the Bifidobacterium shunt, glycolysis, and homolactic fermentation (Spearman rank correlation, BH-corrected P value of ,0.05) (Fig. 7B). Interestingly, these pathways were identified as aging-specific metabolic features by our random forest classifier. Additionally, Gpr41 expression trended toward correlation with fecal acetate (BH-corrected P = 0.085) (Table S4) and propionate (BH-corrected P = 0.058) (Table S4) and was negatively correlated with two taxa at the genus level (Alistipes and Parabacteroides) (BH-corrected P , 0.05) (Table S4). Interestingly, these two taxa were also negatively correlated with fecal SCFA (BH-corrected P , 0.05) (Fig. 6B), indicating a cooccurrence of alterations in the gut microbial community, microbial SCFA, and host GPR41 with aging. Thus, we propose that the gut microbiome shifts as the host ages, resulting in altered microbial carbohydrate metabolism that could in turn impact colon SCFA concentrations, which could feed forward to impact Gpr41/43 expression and host physiology.
As SCFA have been reported to influence multiple facets of inflammation, through their actions both as G-protein-coupled receptor ligands and as inhibitors of histone deacetylases (38), we measured colon gene expression and serum levels of interleukin-6 (IL-6) and tumor necrosis factor alpha (TNF-a) (Fig. 7C and D). A trend of increased Il-6 colon gene expression (P = 0.051) and serum TNF-a levels (P = 0.087) was found in the aged mice, indicating increased inflammation in the aged mice, which is consistent with the general observation of low-grade chronic inflammation status in the elderly (inflammaging) (39).

DISCUSSION
Cooperation between the gut microbiome and the host is finely tuned, with the gut microbiome constantly integrating signals from the host and responding throughout the lifetime of the host (40,41). Host factors such as nutrient metabolism and immune status can impact gut microbiome homeostasis significantly (42). During the aging process, biological and physiological functions systemically decline, leading to a loss of immune and Aging Murine Gut Microbiome Meta-analysis mSystems metabolic homeostasis (3). Thus, we might predict that the gut microbiota would undergo both compositional and functional alterations. Considering the important roles described for the gut microbiome in numerous disease processes, many of which are related to aging (43,44), determining how the composition and function of the gut microbiome change in response to aging, and how these may differ in individuals with healthy versus unhealthy aging, is of great interest. Robust characterization of age-related changes is critical for developing microbiome-based antiaging strategies (e.g., dietary intervention) and identifying therapeutic targets to promote health and longevity in the elderly. In this meta-analysis, we aimed to define aging-related gut microbiome alterations shared across a number of independent mouse studies. Our study found that the aged murine gut microbiome is characterized by increased alpha diversities. Microbial diversity was reported to be decreased with aging in a few human studies (4,5), but this decrease is not consistently observed. Odamaki et al. reported increased diversity in the elderly Japanese population compared to young adults (45), while no difference in alpha diversity was observed between young and elderly Italian populations (11). The variability across human studies might be due to the different patterns of aged gut microbiomes between healthy and relatively less healthy populations as well as betweenpopulation differences in external aging-associated factors impacting the gut microbiome. For example, diet, medication, disease status, and duration of stay in a care facility may all impact the aging microbiome, yet data on these factors are frequently lacking in human aging microbiome research (15,46). Although decreased microbial diversity is associated with several diseases (47), higher diversity does not necessarily imply a healthier microbial community. Additional metrics to describe microbial stability, functional diversity, and components that benefit the host need to be taken into account in aging microbiome research. Thus, the biological meaning of the increased diversity in response to aging requires further investigation.
We further characterized differences in young and old microbial communities. Interestingly, we found a significant age-dependent shift of the microbial community based on unweighted UniFrac distance but not weighted UniFrac distance. As unweighted UniFrac distance takes only absence/presence information into account regardless of the relative abundance (48), the results suggest that rare taxa rather than dominant taxa drove the microbial structure differences between young and old mice. This is consistent with the finding of unaltered F/B ratios between young and aged mice. Instead, we found that a group of rare or subdominant taxa was differentially enriched (e.g., Odoribacter and Mucispirillum) or depleted (e.g., Akkermansia) in response to aging. It is interesting that these taxa are also commonly observed to be altered with aging in human studies (10,11,49,50) and have been implied to be beneficial microbes or pathobionts. Odoribacter has been previously reported as an unfavorable species for metabolism, associated with abnormal lipid and glucose metabolism in the host (51,52). Both Mucispirillum and Akkermansia degrade mucin, yet while Akkermansia is regarded as a beneficial microbe (53), the health relevance of Mucispirillum is still controversial (54). Akkermansia supports gut barrier integrity and contributes to anti-inflammation and antioxidation in aging (55,56). Mucispirillum has been associated with intestinal inflammation in a few studies (57-59), while a protective role of Mucispirillum has been studied in Salmonella colitis in mice (60).
Alterations in rare taxa alone would not necessarily drive significant changes in microbial function. To examine whether aging alters the functional capacity of the gut microbiome, we next examined community-level microbial function using the bioinformatic tool PICRUSt2 to predict microbial metagenome function. We then established a random forest classifier to differentiate microbiomes from young and old animals and to identify aging-specific metabolic features. A random forest classifier is an ensemble of decision trees that uses random selection of different subsets of the features in the training data for constructing the different base classifiers (61). Unlike a single-decision tree, a random forest averages weak learners and avoids putting too much weight on outlier decisions (62). Random forest is thus one of the most robust classification algorithms developed to date and is used extensively in biomedical domains (63). However, when there is a large number of features and the proportion of truly informative features is small, the performance of a random forest classifier is likely to decline due to the low probabilities of choosing truly informative features by random sampling. Therefore, we applied an iterative random forest algorithm to construct the random forest classifier (64). The algorithm starts with the entire set of metabolic features in the training data set and stores the importance (mean decrease in Gini impurity) of the metabolic features. In the next iteration, a new weighted random forest is grown with weights set to the feature importance from the previous iteration. Ten iterations were performed in our study. Considering the limitations of using predicted metagenomic function rather than direct metagenome sequencing to construct the random forest classifiers, we included one 16S rRNA gene sequencing and one whole-genome sequencing data set generated by our laboratory as external data sets to validate the final random forest model. The performance of our random forest classifier in differentiating between microbiomes from young and old mice was good, with obtained AUC values of between 0.75 and 0.97 for the test data set and the two external data sets. Thus, the features selected by the random forest classifier appear to be informative and reliable aging-specific metabolic features.
Interestingly, we found that 50% of the 20 most important features identified by the random forest classifier are involved in carbohydrate metabolism. Bacterial carbohydrate metabolism is a basic biological process providing energy and precursors for biosynthetic pathways. Additionally, metabolites derived from microbial carbohydrate metabolism, primarily short-chain fatty acids (SCFA), exert biological effects on the host (38). The altered carbohydrate metabolism potential that we observed in aged mice potentially affects the production of SCFA. In fact, we found that fecal acetate, propionate, and butyrate were significantly decreased in the aged mice, which is consistent with data from other published mouse studies (26,27). GPR41 and GPR43 are the primary receptors that sense SCFA and mediate the biological effect of SCFA on host health (37). GPR41 prefers the longer SCFA (butyrate $ propionate . acetate), while GPR43 responds more strongly to the shorter SCFA (acetate $ propionate . butyrate) (37). The decreased Gpr41 expression and increased Gpr43 expression in aged mice suggest that host-microbiome communication could be altered with aging both by changes in the abundance of SCFA and by the decreased ability to sense longer SCFA (e.g., butyrate and propionate). The correlation of colonic Gpr41 expression in the host with the abundance of gut microbes, microbial carbohydrate degradation and fermentation metagenomic potential, and fecal SCFA suggests that changes in the gut microbiome-SCFA with aging may in turn impact GPR41 and downstream host pathways. One potential downstream effect of altered SCFA-GPR41/43 signaling is inflammation, and thus, we examined two cytokines reported to be increased with aging. Despite a relatively small data set, we detected a trend toward increased Il-6 gene expression in the colon as well as increased serum TNF-a in aged mice, suggesting an increased inflammatory state in aged mice. However, the extent to which changes in the gut microbiome-SCFA-GPR41 axis with aging impact host physiology and health remains to be elucidated.
There are several limitations of our study. To generate comparable data sets, a closed-reference OTU-picking process was chosen. This method restricts OTUs identified within the reference database (Greengenes V13-8), which excludes OTUs that are not present in the database and biases toward lower microbial diversity. Moreover, taxonomic assignments between studies may be influenced by the differences in sequencing lengths and regions between studies, which potentially introduce noise across the data sets (65). In addition to these technical issues, host-related covariates could also contribute to variations across the studies, which we were not able to control due to the missing individual-level metadata of aging-associated biological factors (e.g. metabolic and inflammation status [16,66]). Furthermore, aging is a gradual process, and gut microbiome alteration with aging is likely to also be a gradual process in mice. However, due to the limited data available from published studies, we were able to include samples at only two ages bookending the spectrum of aging (young versus old). Additionally, we defined a relatively wide range for old (18 to 26 months of age). Thus, we lacked the ability to detect differences between moderately and extremely old mice and potentially introduced variability in our aged cohort that could decrease our ability to detect true differences between young and old mice. Finally, we cannot exclude the possibility that the observed correlation of host Gpr41 expression with the gut microbiome is driven by age-related changes in both. Due to the limited age variation of our samples, we were unable to effectively control for age as a covariate in our Spearman correlation analysis. In the future, we hope to expand this work to include a variety of ages across the life span, which would facilitate investigations of whether correlations between changes in gut microbiome features and alterations in host physiology are independently correlated.
Despite these limitations, our study revealed both compositional and functional alterations in the gut microbiome with aging. Altered carbohydrate metabolism was a key feature of aging-related changes in the gut microbiome and was accompanied by decreased fecal SCFA. On the host side, colonic Gpr41 expression was significantly correlated with microbial carbohydrate degradation and fermentation pathways. These findings show that compositional changes in the gut microbiome with aging correlate with functional changes, leading to measurable effects on key metabolites, as well as age-related changes in host receptors for those metabolites and imply cross talk between the gut microbiome, microbial SCFA, and host GPR41-dependent pathways in response to aging. We hope that the results from this meta-analysis across 5 independent studies will be a foundation for future work to further advance our understanding of the aging microbiome and host interaction from a putative correlation to a causative relationship and, ultimately, translation to clinical applications such as microbial manipulation to improve healthy aging and prevent age-related diseases.

MATERIALS AND METHODS
Data sets. Literature searches of PubMed and Google Scholar were executed in July 2020 using the following keyword terms: "aging microbiome" OR "aging microbiota" OR "aged microbiome" OR "aged microbiota" OR "old microbiome" OR "old microbiota" AND "murine" OR "mice" OR "mouse." The search yielded 35 potential studies for inclusion. We excluded studies where raw sequencing data were not publicly available, the sample number was ,3 per age group, the average sequencing depth was ,10,000 reads per sample, or the study groups did not meet our definition of young (1 to 6 months old) or old (18 to 26 months old). We also excluded studies of genetically modified or diseased animals (see Fig. S1 in the supplemental material).
After filtering, 5 studies that fulfilled our criteria were included in the meta-analysis. One study included data from both female and male C57BL/6 mice (18), and these were analyzed separately for a total of 6 data sets: 3 from female (16-18) and 3 from male (18-20) C57BL/6 mice. Within each data set, we further filtered out individual samples from animals that had a low sequencing depth (,2,000 reads) or were manipulated in any way (i.e., dietary intervention). The raw sequencing data were downloaded from the SRA and ENA databases by accession number, as listed in Table S1. As external data sets for validation of the random forest classifier, we included one fecal whole-genome metagenomics data set from male C57BL/6 mice sequenced as described below (NCBI SRA accession number PRJNA739153) and one 16S rRNA gene data set from male CB6F1 mice previously sequenced by our laboratory (You X, Yan J, Herzog J, Campbell R, Hoke A, Hammamieh R, Sartor RB, Kacena MA, Chakraborty N, Charles JF, manuscript in preparation) (NCBI SRA accession number PRJNA737742).
Thus, in total, 3 female C57BL/6 data sets (n = 32) and 3 male C57BL/6 data sets (n = 102) were included in the meta-analysis and the generation of the random forest classifier. External data sets for validation of the random forest classifier included one male CB6F1 16S rRNA gene sequencing data set (n = 16) and 1 male C57BL/6 whole-genome metagenomics data set (n = 20).
OTU picking. Due to the variant amplicons of the 16S region that were sequenced across studies, we chose a closed-reference OTU-picking process to ensure concordance between studies. The raw sequencing data were imported into QIIME2/2018.11 for data processing (67). DADA2 was used to denoise reads (68). The reads were trimmed such that the mean quality score was ,30 for four consecutive reads or the minimal sequencing length. For paired reads, reads were merged if the overlapped reads were .20 bp after truncation. Otherwise, only forward reads were used for the following analysis. Due to the potentially insufficient removal of the adapter/index/primer in the raw sequencing data, which would interfere with DADA2's error model, we also trimmed 16 nucleotides from the 59 end of each read. The quality-controlled, dereplicated reads were then subjected to OTU picking by q2-vsearch against the Greengenes V13-8 database at 97% identity.
Diversity analysis. Six OTU tables generated from the q2-vsearch clustering process were combined. OTUs with fewer than 5 total reads or that presented in only one sample were filtered out. Samples were rarefied to 10,000 reads by the rrarefy function in the R package Vegan (2.5-6) (69). Shannon, Simpson, and Chao1 indices were calculated by the diversity and estimateR functions in Vegan. Faith's phylogenetic distance was calculated by the R package Picante (1.8.2) (70). Welch's t test was performed per study analysis using the function t.test (stats, R package). The combined studies were analyzed using a linear mixed-effects model by the lmer function in the R package LmerTest (3.1-3) (71) using the formula log 2 FC ; age 1 (1jstudy). The same formula was also applied to the enrichment analysis of centered log-ratio normalization of OTUs and predicted PICRUSt2 metagenomic functional pathways (PICRUSt2 method [see below]). Weighted UniFrac and unweighted UniFrac distances were calculated by the R package Phyloseq (1.32.0) (72). PERMANOVA with 999 permutations (vegan::adonis) was performed for the differences in beta diversities. For combined studies, PERMANOVA was performed with the option strata to control for study.

Aging Murine Gut Microbiome Meta-analysis mSystems
Fecal DNA extraction and whole-genome microbial sequencing. Individual fecal samples were collected from 10 young (2-month-old) and 10 old (26-month-old) C57BL/6 male mice (NIA Aged Rodent Colony, USA), under IACUC approval 20110 (Indiana University). Genomic DNA was isolated from fecal samples using a QIAamp PowerFecal Pro DNA kit (Qiagen, MD, US) according to the manufacturer's protocol. Fecal DNA was diluted to 10 ng/mL and submitted to the Broad Institute (MA, US) for high-output whole-genome microbial 150-bp paired-end sequencing. After sample DNA quality assessment, sequencing libraries were prepared using the Illumina Nextera XT DNA library prep kit. Two positive controls and one negative control were included in the sequencing run. In total, 1126.04 million reads (56.30 million 6 12.42 million reads per sample) were generated (Fig. S4).
Whole-genome OTU picking and metagenomic functional analysis. Unmapped BAM files were received from the Broad Institute Genomics Center. Picard (2.6.0) was used to convert BAM files to fastq files. Sequencing quality was checked by FastQC (0.11.8). Raw reads were adapter trimmed by Trimmomatic (0.39) (73) and decontaminated of the host genome (mouse GRCm38.p6) mapped by Bowtie2 (2.4.2) (74) as part of the Kneadata (0.7.4) pipeline using the default parameters (https://github .com/biobakery/kneaddata). To consistently apply the analysis pipeline to 16S data sets, we used an OTU-picking process to generate the OTU table. Filtered reads were selected by alignment against the Greengenes V13-8 database via SortMeRNA (v.4.3.4) (75). OTU picking on the resulting reads was then performed by q2-vsearch against the Greengenes V13-8 database at 97% identity. For metagenomic functional analysis, filtered reads were processed by the HUMAnN3 (3.0.0.alpha.3) pipeline for the metagenomic functional analysis, including the following steps: (i) taxonomic profiling by MetaPhlAn3 (3.0.7) (76), (ii) indexing of the identified species pangenomes in the ChocoPhlAn database and alignment of reads to the reference pangenomes by Bowtie2, (iii) translated searching of unmapped reads against the UniRef90 protein reference database using DIAMOND (2.0.6) (77), and (iv) mapping of the resulting gene lists to the MetaCyc pathways. The normalized relative abundances of OTUs and MetaCyc pathways were then used as an external data set for the validation of random forest classifiers as described below. Moreover, the taxon table obtained by MetaPhlAn3, which has a better taxon resolution than OTU picking, was used to perform Spearman correlation analysis.
Random forest classifier. PICRUSt2 (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States) (78) was used to predict metagenomic functions based on the six 16S rRNA gene data sets described above in "Data sets." The data sets of MetaCyc pathways and OTUs were used as predictor variables for the random forest classifier. The data set was randomly divided into training (n = 89) and testing (n = 45) data sets. The iterative random forest (iRF) algorithm was performed using the R package iRF (64). Ten iterations were performed, of which the random forest with the smallest out-ofbag (OOB) error rate was chosen.
The receiver operating characteristic (ROC) curve and the area under the ROC curve (AUC) were calculated for the testing data set by the prediction and performance functions in the R package ROCR (79). We then applied the same functions to the two external data sets generated in our laboratory and as described above in "Data sets" to validate the performance of the random forest classifier.
Fecal samples were homogenized in 400 mL of water for 2 min on ice and centrifuged, and the supernatant was collected. For derivatization purposes, 20 mL each of AABD-SH, TPP, and DPDS were added to the supernatant and vortexed at room temperature for 5 min. The reaction mixture was dried under a vacuum and reconstituted with 150 mL of methanol containing lactic acid-13 C 3 as an internal standard. The samples were centrifuged, and the supernatant was diluted 10-fold.
Three microliters of the prepared sample was injected onto a Kinetex 2.6-mm, 100-Å, 100-by 3.0-mm polar C 18 column (Phenomenex, CA, USA) using a SIL-30 AC autosampler (Shimadzu) connected to a high-flow LC-30AD solvent delivery unit (Shimadzu) and a CBM-20A communication bus module (Shimadzu) online with a Qtrap 5500 instrument (Sciex, MA, USA) operating in positive-ion mode. A binary solvent comprising water with 0.1% formic acid (solvent A) and acetonitrile with 0.1% formic acid (solvent B) was used. The extracted metabolites were resolved at a 0.4-mL/min flow rate. MS was conducted in positive-ion mode with a turbo ion spray voltage of 5,500 V, using 20 lb/in 2 of curtain gas, 50 lb/in 2 of nebulizer gas, and 50 lb/in 2 of drying gas at a temperature of 400°C. LC separation was performed using mobile phase A (0.1% formic acid in water) and mobile phase B (0.1% formic acid in acetonitrile) at a flow rate of 400 mL/min and a temperature of 40°C. The separation gradient was as follows: 0% B at 0 min, 0 to 100% B in 8 min, 100 to 0% B in 1 min, and 0% B in 1 min. A collision energy of 15 V was used for multiple-reaction monitoring (MRM), and LC-tandem MS (LC-MS/MS) data were analyzed by Analyst 1.5.2 software (AB Sciex). The peak area of the lactic acid-13 C 3 isotope-labeled internal standard was used to normalize SCFA and processed using MultiQuant 3.0.3 software. The column was conditioned using the pooled quality control (QC) samples initially and after every 6 sample injections to monitor shifts in signal intensities and retention times. A blank solvent was run before and after the pooled QC samples to minimize carryover effects. The concentration data were calculated from standard curves performed by QC-based locally estimated scatter plot smoothing (LOESS) signal correction (QC-RLSC) followed by normalization to the input fecal weight.
Colon RNA extraction, cDNA synthesis, and qRT-PCR. Colon tissues were snap-frozen in liquid nitrogen and stored at 280°C until analysis. Total RNA was isolated by tissue homogenization with a bullet blender bead beater (Navy Rino RNA lysis kit; Next Advance) followed by an RNeasy minikit (Qiagen) according to the manufacturer's protocol. DNA contamination was removed by DNase I treatment (Thermo Scientific). cDNA synthesis was performed using a high-capacity cDNA reverse transcription kit (Applied Biosystems). Quantitative real-time PCR (qRT-PCR) was performed using Fast SYBR green master mix (Applied Biosystems). The primer sequences of the target genes are listed in Table S2. The reactions were performed on a StepOne Plus real-time PCR machine (Applied Biosystems, USA). Results were normalized to the housekeeping gene hypoxanthine phosphoribosyltransferase (HPRT) and calculated by the 2 2DDCT method.
Serum TNF-a and IL-6 analysis. Blood samples were collected by cardiac puncture. Serum samples were separated by clotting the blood samples for 30 min, followed by spinning for 10 min at 12,000 Â g. Serum TNF-a and IL-6 were measured using a murine TNF-a standard 2,29-azinobis(3-ethylbenzthiazolinesulfonic acid) (ABTS) enzyme-linked immunosorbent assay (ELISA) development kit and a murine IL-6 standard ABTS ELISA development kit (Peprotech) according to the manufacturer's protocol.
Statistical analysis. Unless otherwise specified, statistical analysis was performed with GraphPad Prism 8. Data were represented as means 6 standard deviations (SD). An unpaired t test (parametric data) or a Mann-Whitney U test (nonparametric data) was performed for fecal SCFA, colon gene expression, and serum TNF-a and IL-6. Spearman rank correlation was performed between colon gene expression and gut microbial carbon degradation/fermentation pathways, between colon gene expression and fecal SCFA, and between taxa and fecal SCFA. The adjusted P value was determined by the Benjamini-Hochberg method using the function p.adjust in the base R package for P values obtained from the Spearman rank correlation and enrichment analyses. A P value of ,0.05 was regarded as statistically significant.
Data availability. Raw sequencing data for each study can be accessed as listed in Table S1. The sequencing data generated in our laboratory has been deposited in the NCBI SRA under accession numbers PRJNA737742 (16S rRNA gene) and PRJNA739153 (whole-genome metagenomics).

SUPPLEMENTAL MATERIAL
Supplemental material is available online only.


Methionine regulates self-renewal, pluripotency, and cell death of GIC through cholesterol—rRNA axis

Background Glioma-initiating cells (GICs) are the source of glioma cells that can self-renew, have pluripotency, and are treatment-resistant, so are the starting point for relapse and eventual death despite multimodality therapy. L-[methyl-11C] methionine PET has observed high accumulation at the time of recurrence, it is important to understand the mechanism of tumor cell activation caused by the reorganization of methionine metabolism. Methods We cultured cells in methionine-deprived culture medium for comprehensive analysis. Based on the obtained results, the possible target molecules were chemically inhibited and the respective markers were analyzed. Results Methionine depletion markedly decreased proliferation and increased cell death of GICs. Decreased S-adenosyl-methionine, which is synthesized intracellularly by catalyzed by methionine adenosyltransferase using methionine, triggered the following: (i) global DNA demethylation, (ii) hyper-methylation of signaling pathways regulating pluripotency of stem cells, (iii) decreased expression of the core-genes and pluripotent markers of stem cells including FOXM1, SOX2, SOX4, PROM1, and OLIG2, (iv) decreased cholesterol synthesis and increased excretion mainly through decreased SREBF2, and (v) down-regulation of the large subunit of ribosomal protein configured 28S and ACA43, small nucleolar RNA guiding the pseudouridylation of 28S rRNA, which is essential for translation. In addition, inhibition of cholesterol synthesis with statin resulted in a phenotype similar to that of methionine depletion and decreases in stem cell markers and small nucleolar RNA ACA43. Moreover, suppression of FOXM1 decreased stem cell markers such as SOX4 and PROM1. The gene expression profile for cholesterol production was obtained from the Ivy Glioblastoma Atlas Project database and compared between tumor cells with relatively low methionine levels in areas of pseudopalisading arrangement around necrosis and tumor cells in the infiltrating region, showing that cells in the infiltrating region have higher capacity to produce cholesterol. Conclusions Methionine metabolism is closely related with self-renewal, pluripotency, and cell death in GICs through modification of cholesterol biosynthesis, especially in the SREBF2-FOXM1 and ACA43 axis with modification of rRNA. Supplementary Information The online version contains supplementary material available at 10.1186/s12885-022-10280-5.

Introduction
Glioblastoma multiforme is the most common and deadly type of primary brain tumor with a median overall survival of less than 15 months despite aggressive treatment [1]. Glioblastoma multiforme can produce glioma Open Access *Correspondence: akatoyik@med.miyazaki-u.ac.jp cells and consequently self-renew, so glioma-initiating cells (GICs) are considered to be the root stem cell for tumorigenesis and recurrence. Cancer stem cells including GICs retain "plasticity" that changes dynamically as a result of metabolic reorganization.
Malignant transformation or recurrence of glioblastoma multiforme can be detected as increased methionine accumulation on PET before contrast enhancement is seen on MRI. Therefore, detailed understanding of the methionine metabolism of GICs may be important for tumor control. Methionine concentrations are expected to be high in areas of sparse tumor cells, but relatively low in areas of dense lesions, such as in the peri-necrotic pseudopalisading region, due to low blood supply. Therefore, glioblastoma multiforme may develop with nonuniform concentration of methionine within the same lesion. We hypothesized that reorganization of the methionine metabolism in the presence of restricted common source of methyl donors could modify the character of tumor cells as reflected by epigenetic status, and gene and protein expression profiles.
The present study showed that the methionine metabolism in GICs was closely related to cholesterol metabolism and rRNA expression, resulting in change in pluripotency, progression, and avoidance of cell death.

Cell lines and culture conditions
Human GICs (MZGC1, MZGC2, MZGC3, MZGC4, MZGC5, MZGC6, MZGC7 and MZGC8 cells) were isolated and established as cell lines from the surgical tissues of patients treated in our institution from March 2015 to April 2021. GICs were cultured in Dulbecco's modified Eagle's medium/nutrient mixture F12 Ham's liquid supplemented with SUBSER-ESrP ® (Cell Science & Technology Institute, Inc., Sendai, Japan), CELRENA ® medium (Cell Science & Technology Institute, Inc.), in 2-hydroxyethyl methacrylate -coated T75 flasks in a humidified incubator at 37 °C under an atmosphere of 5% CO 2 and 95% air. For the methionine depletion studies, cells were plated on laminin-coated 10-cm culture dishes and cultured with complete CELRENA ® medium or methionine deficient CELRENA ® medium (custom medium CST1 MFG No. cn2243a0) for 48-72 h. Cells were dispersed onto laminin-coated 96-well plates, replaced with treated media on day 1, and time-lapse photography was performed with Incucyte ® (Sartorius, Osaka, Japan). Simvastatin and FDI-6 were purchased from Selleck (Tokyo, Japan).

Dye exclusion test
Dye exclusion testing was performed to analyze cell viability. Fixed amounts of cells were seeded and cultured in medium with (control) or without methionine for 3, 5 and 7 days. Harvested cells were suspended in phosphate buffered saline-containing trypan blue and then examined to determine the percentage of cells with blue clear cytoplasm (nonviable cells) versus total cells using a cell counter (LUNA, Logos Biosystems, Annandale, VA, USA). The relative rate against the control sample was calculated.

Cell cycle analysis and detection of surface marker of GICs
Cells with or without treatment of methionine depletion/statin addition were collected using Accutase ® solution (Sigma, USA). The Cell Cycle Phase Determination Kit (Cayman Chemical, Ann Arbor, MI, USA) was used for cell cycle analysis. In brief, collected cells were rinsed twice with buffer, then fixed at − 20 °C overnight. Cells were washed twice with ice-cold phosphate buffered saline, and stained with propidium iodide/RNase staining buffer solution in the dark for 30 min at room temperature. Then, cells were analyzed with a flow cytometer (Guava ® EasyCyte ™ Mini, Luminex Japan, Tokyo, Japan). A histogram of the cell cycle distribution was generated from 5000 events per sample and data were analyzed using Guava ® Cell Cycle software. To detect CD133, the putative marker of GICs, by flow cytometry (FCM), 1 × 10 6 cells derived from tumorspheres were analyzed using fluorescent (FITC) labelled antibodies (CD133/1 (AC133)-VioBright FITC, Miltenyi Biotec, Germany) and BD FACSCalibur (BD Biosciences, USA).
RNA extraction including microRNA (miRNA), cDNA synthesis, and quantitative PCR (qPCR) analysis RNA from cultured tumor cells was extracted with the mirVana miRNA Isolation Kit (Ambion, Thermo Fisher Scientific K.K., Tokyo, Japan) or RNeasy Plus Mini Kit (QIAGEN, Germantown, MD, USA). For cDNA synthesis, RNA was reverse transcribed from random hexamers using the SuperScript ™ VILO ™ cDNA Synthesis Kit (Invitrogen, Life Technologies, Thermo Fisher Scientific, Waltham, MA, USA). Real-time qPCR was then performed in triplicate on the StepOne Plus or Quant3 (Applied Biosystems, Thermo Fisher Scientific) using SYBR ™ Green Realtime PCR Master Mix (Applied Biosystems, Thermo Fisher Scientific) to determine the mRNA levels. PCR was performed using a 20 μl volume containing 2 µl cDNA, 300 µM of each primer, and 10 µl of 2 × PCR master mix under the following conditions: 95 °C for 10 min followed by 40 cycles of 95 °C for 15 s and annealing/extension at 60 °C for 1 min. The data were normalized to the amount of human 18S rRNA, and the values are represented as the mean ± SD of 2 −ΔΔCt in a triplicate assay.
The primers used in this study were described in Supplementary Table S2.

Colony formation assay
Cells were seeded into laminin-coated 12-well plates in triplicate at a density of 500 cells/well in 2 ml of medium with or without treatment medium. After 120 h, formed colonies were trypsinized, reseeded and cultured for further 72 h. After 192 h, the cell colonies were stained for 15 min with a solution containing 0.5% crystal violet and 25% methanol, followed by three rinses with tap water to remove excess dye. Each well was observed with Olympus IX71 microscope (Olympus, Tokyo, Japan).

Gene and miRNA expression analysis
Gene and miRNA expression was analyzed with a Gene-Chip ™ System with a Human Genome Clariom ™ D Array and GeneChip ™ miRNA 4.0 Array (Thermo Fisher Scientific) according to the manufacturer's instructions. mRNA expression analysis used amplification and biotin labeling of fragmented cDNA with the biotin labeling system (GeneChip ™ WT PLUS Reagent Kit, Thermo Fisher Scientific) in duplicate. miRNA expression analysis used biotin-labeled total RNA 1 µg including miRNA from tissue using the FlashTag ™ Biotin HSR RNA Labeling Kit (Affymetrix ® , Thermo Fisher Scientific). Labeled probes were hybridized to the Human Genome Clariom ™ D Array (Thermo Fisher Scientific) or GeneChip ™ miRNA 4.0 Array (Thermo Fisher Scientific) with the GeneChip ™ Scanner 3000 7G (Affymetrix ® , Thermo Fisher Scientific). Expression data of mRNA were extracted from image files using GeneSpring GX 14.9.1 (Agilent Technologies, Santa Clara, CA, USA). Normalization and expression value calculations were performed using the GeneSpring system. Signal values of miRNA were calculated using the Transcriptome Analysis Console software (Affymetrix ® , Thermo Fisher Scientific). Statistical analysis was performed using two-tailed, unpaired t-tests. Any change ≥ twofold with P < 0.05 was considered statistically significant.
Enrichment analysis used Metascape (http:// metas cape. org/) and Gene Set Enrichment Analysis (GSEA) (http:// www. broad. mit. edu/ GSEA). Significance of biological processes was determined through P-values calculated on a hypergeometric distribution (log10). Metascape was used to conduct meta-analysis with generation of heat maps of gene ontology terms that hierarchically cluster together under experimental conditions. Statistical criteria for a differentially expressed gene included a change greater than 1.5 fold, and an overall false discovery rate smaller than 5% (change ≥ ± 1.5 fold and false discovery rate ≤ 0.05) with Benjamini-Hochberg correction (q) to account for multiple comparisons. GSEA using MSigDB (v7.4) was also used to determine the significance of a pre-defined gene set by comparing the correlation between expression and class distinction with other random situations. The significance threshold was set at nominal P value < 0.05 or false discovery rate q value < 0.25.

Measurement of cholesterol, S-adenosyl-methionine (SAM), and S-adenosylhomocysteine (SAH) in tissue using ELISA
Measurement of cholesterol, SAM, and SAH used glioblastoma stem cells cultured with CELRENA medium or methionine-deprived CELRENA medium for 72 h and 10 × 10 6 cells/cell pellets were homogenized by sonication in ice-cold phosphate buffered saline followed by centrifugation at 10,000 g for 15 min. Intracellular cholesterol, SAM, and SAH concentrations were measured using the Total Cholesterol Assay Kit (Colorimetric, Cell Biolabs, Inc., San Diego, CA, USA) and SAM and SAH ELISA Combo Kit (Cell Biolabs, Inc.) in accordance with the manufacturer's instructions under the indicated conditions. The standards were generated using the cholesterol standard and SAH/SAM-bovine serum albumin conjugate supplied in the kit. After the reaction was stopped, the OD450 and OD620 values were read using a SpectraMax iD3 Multi-Mode Microplate Reader (Molecular Devices, San Jose, CA, USA). SAH/SAM ratios were calculated as the ratio of SAM to average SAH using the intracellular concentration (µg/mL) to determine the relative methylation potential.

Reduced Representation Bisulfite Sequencing (RRBS) and data analysis
RRBS libraries with single MspI digestion were constructed for glioblastoma stem cells. Briefly, 1 ug of genomic DNA was digested with 20U of MspI enzymes in 18 ul reaction mixtures at 37 °C for 2-3 h. After purification, the digested products were blunt-ended, and then dA added, followed by methylated-adapter ligation using a NEXTflex Bisulfite-Seq Ligation Kit. To obtain DNA fractions of > 100 bp, the MspI-digested products were purified using Agencourt AMPure XP (Beckman Coulter). Bisulfite conversion was conducted using an EX DNA Methylation Gold Kit following the manufacturer's instructions. The final libraries were generated by PCR amplification using NEXTflex Bisulfite-Seq Ligation Kit and sequenced using an Illumina HiSeqX to generate 150 bp pair-end reads (Supplementary Table S1).
CpG sites in the resulting RRBS data were then interrogated for methylation patterns and differential methylation (q value < 0.05 and methylation percentage difference of at least 25%) using the methylKit. The differential methylation data were then queried for differentially methylated regions (DMRs) using annotatr (https:// bioco nduct or. org/ packa ges/ relea se/ bioc/ html/ annot atr. html) and GeneDMRs, R package for Gene-based DMR analysis (https:// www. biorx iv. org/ conte nt/ 10. 1101/ 2020. 04. 11. 03716 8v1. full). The resulting DMR outputs were also visualized in the genome browser or in the Integrative Genomics Viewer (IGV) program (https:// softw are. broad insti tute. org/ softw are/ igv/). To establish the relationship between histone mark, insulators, transcription factors, and DMRs, the RRBS data conjunction with the data acquired from the Database of ChIP-Atlas Peak browser (https:// chip-atlas. org/ peak_ brows er) were projected on to the IGV.

In silico analysis
Gene expression data obtained by the RNA sequencing technique involved in cholesterol metabolism in different anatomical regions (infiltrating tumor, pseudopalisading cells) of the glioblastoma multiform were identified using the Ivy Glioblastoma Atlas Project (Ivy GAP) database (http:// gliob lasto ma. allen insti tute. org/).

Statistical analysis
Statistical analyses were performed using GraphPad Prism 7 (GraphPad Software, La Jolla, CA, USA). Data are expressed as means ± SD of experiments performed in triplicate. One-way analysis of variance with Tukey's post-hoc test or unpaired Student t-test was performed for comparisons of the quantitative data between groups. P < 0.05 was considered to indicate a statistically significant difference.

Methionine depletion decreased cell proliferation and increased cell death
Clonogenic assay showed marked decreases in cell growth and sphere formation, suggesting reduced selfrenewal capacity induced by the depletion of methionine (Fig. 1A). Fluorescence-activated cell sorter analysis Quantitative analysis of colony areas in response to methionine depletion. B Fluorescence-activated cell sorter analysis to detect cell cycle distribution. Propidium iodide (PI) staining assay using a Cell Cycle Phase Determination Kit was performed on MZGC1 cells. Cell cycle analysis showed that methionine depletion for 3 days increased sub-G1 and decreased S, G2M populations. C Dye exclusion test showed that methionine depletion for 3 days increased dead cell rate. D Fluorescence-activated cell sorter analysis using annexin-V kit showed that methionine depletion induced apoptotic cell death and dye exclusion test were also performed to assess cell death or cell cycle arrest. Methionine depletion for 3 days resulted in increased sub-G1 population and cell death ( Fig. 1B and C). Further fluorescence-activated cell sorter analysis using annexin-V showed apoptosis induced by methionine depletion, although the time to apoptosis varied from cell to cell (Fig. 1D, Supplementary Figure  S1).

Methionine depletion impaired SAM-dependent methyltransferase pathway
A comprehensive gene expression analysis was performed using a microarray to detect the characteristics of changes in specific molecular networks and canonical signaling pathways induced by modified methionine metabolism. A total of 213 genes were upregulated in the methionine withdrawal group, which were enriched in 17 Metascape terms/pathways, including 1) SAMdependent methyltransferase pathway, 2) hallmark TNFA signaling via the nuclear factor-kappa B (NFKB), 3) regulation of cytokine biosynthetic process, 4) neutral amino acid transmembrane transporter activity, 5) circadian rhythm, 6) ribonucleoprotein complex biogenesis, 7) negative regulation of protein modification process, and 8) positive regulation of cholesterol efflux ( Fig. 2A, Supplementary Figure S2). In contrast, 115 genes downregulated in the methionine withdrawal group were enriched in terms, regulation of cholesterol metabolic process and neural crest cell migration. SAM is a methionine derivative that is essential for the catalytic reaction of methylase. Intracellular SAM levels were reduced by methionine depletion (Fig. 2B). SAM is synthesized intracellularly by methionine adenosyltransferase (MAT) using the essential amino acids methionine and ATP as substrates. SAM is a methyl group donor and the dependent pathways include MAT2A, MAT2B, DNMT1, DNMT3A, and DNMT3B. qPCR and expression array were performed to examine whether methionine depletion in GICs showing that methionine depletion significantly activated hallmark TNFA signaling via the NFKB pathway. GSEA of differential expression between control and methionine depletion states demonstrated that methionine-depleted GICs were enriched in hallmark TNFA signaling via the NFKB pathway. Heat map of top 20 genes from GSEA analysis is shown on the right. Red column corresponds to upregulated genes and blue column corresponds to downregulated genes could functionally regulate histone and DNA methylase. Methionine depletion up-regulated MAT2A and MAT2B in time-dependent manners (Fig. 2C). DNMT3B, DNMT3A, TET1, and TET2 were down-regulated (Fig. 2D). As shown above, hallmark TNFA signaling via the NFKB is activated in the methionine-depleted state, as we confirmed in GSEA (Fig. 2E). Cell death via apoptosis might result from the activated TNFA signaling pathway (Figs. 1B-D, 2A, and E).

Methionine depletion leads to pathway-dependent changes in methylation status
We tried to assess the relationship between transcriptome and DNA methylation status. Epigenetic analysis using RRBS revealed that depletion of methionine promoted genome wide hypomethylation in both promoter and gene body in many cases (Fig. 3A). We analyzed DMRs to identify DMR networks associated with affected pathways using the R package "GeneDMRs. " GeneDMRs can separately analyze the methylation levels within the gene, gene body, CpG island regions (CpGi) mainly localized in the promoter, and overlapping regions (Fig. 3B). Although there is some degree of variation between cell lines (Fig. 3A), analysis of raw data using GeneDMRs allowed statistically significant pathways to be picked up according to the degree of methylation. Pathway analysis of CpGi revealed that signaling pathways such as regulating pluripotency of stem cells, porphyrin and chlorophyll metabolism, retinol metabolism, regulation of actin cytoskeleton, and MAPK signaling pathway are hypermethylated (Fig. 3C) in the presence of methionine depletion. In contrast, regulation of lipolysis in adipocytes, circadian entrainment, neuroactive ligandreceptor interactions, glutamatergic synapse, cholinergic synapse, and phosphatidylinositol signaling system pathways were hypomethylated. CpGi analysis often classifies pathways as either hyper-or hypomethylated, but often overlap by gene body analysis (Supplementary Figure  S3). Pathways classified as either by gene body analysis included focal adhesion and Kyoto Encyclopedia of Genes and Genomes glioma pathway which were hypomethylated (Supplementary Figure S3).

Methionine depletion changed GIC pluripotency
RRBS demonstrated the regulating pluripotency of stem cell pathways (Fig. 3C), so the expression of genes related to GIC pluripotency was investigated. Expression array analysis observed down-regulation of OLIG2 gene expression, one of the four core genes of GICs [3,4], and down-regulation of FOXM1, PROM1, and other marker/maintenance genes of GICs (Fig. 4A). qPCR, FCM and immunoblotting experiments confirmed down-regulation of these markers ( Fig. 4B and C). We investigated epigenetic changes of those genes. The CpGi in the promoter and genebody region of PROM1 gene are mainly hypomethylated by the removal of Methionin, while the shore region is mainly hypermethylated. Interestingly, almost the same locus of the gene body was changed in three GICs (Fig. 4D).
Next, we compared the expression levels of genes of cholesterol metabolism within the tumor mass (perinecrotic pseudopalisading area) with the infiltrating area using Ivy GAP database. The genes related to de . FOXM1 binds to intron of PROM1 gene body, whose histone mark is H3K27ac. Hypermethylated (pink) and hypomethylated (green) in the methionine depleted compared to control states are indicated. Hg19 used as reference genome novo synthesis of cholesterol were predominantly upregulated in infiltrating tumor cells compared with tumor cells in the perinecrotic pseudopalisading area (Fig. 5E, Supplementary Figure S4). However, ABCA1, which acts as the efflux pump of cholesterol, showed the opposite distribution. Simvastatin was added to assess whether chemical inhibition of HMGCR had same effect with methionine depletion. Simvastatin decreased colony formation, showing a very similar phenotypic change to methionine depletion (Fig. 5F,  Figure S6A, B and S10C). These results indicated that methionine metabolism is deeply involved in the survival, proliferation, and pluripotency of GICs through cholesterol metabolism.
IGV investigation of the methylation status of each gene found no significant change in gene body or CpGi as in HMGCS1 (Fig. S7A). We think that the cholesterol metabolic genes might be mainly regulated by another mechanism. We analyzed the transcription factor candidates and insulator binding sites that are likely to bind to the promoter region from the ChIPseq database. CpGi in the promoter of HMGCS1 had no DMRs, but CTCF and SREBF2 binding are reported (Fig. S7A). Since SREBF2 is known to be an important regulator of cholesterol metabolism, we checked the expression of SREBF2. SREBF2 was down-regulated by methionine depletion (Fig. 5G). Moreover, we found multiple genes in the Horton SREBF target gene set were down-regulated by methionine depletion (Fig. 5H). Among the top 10 affected genes related to cholesterol biosynthesis (Fig. 5C), HMGCS1, FDPS, IDI1, SC5D, TM7SF2, and HSD17B7 are SREBF target genes. Therefore, we thought that SREBF2 was one of the key molecules linking methionine and cholesterol metabolism. To assess the methylation status of SREBF2, we applied the data of RRBS to IGV. DMRs were mainly found in the gene body including 3'-untranslated region (Fig. S7B).
We then considered whether cholesterol metabolism is related to the pluripotency of GICs. FOXM1, a maintenance gene of GICs, is also associated with cholesterol metabolism [6]. Thus, we investigated the relationship between cholesterol metabolism and FOXM1. FOXM1 was down-regulated with SREBF2 by methionine depletion (Figs. 4A, B, and 5G). SREBF2 increased transiently in the short time after statin administration with feedback that cholesterol production was reduced, but then decreased in the longer term (Fig. 5I). FOXM1 was up-regulated together with upregulated SREBF2 (Fig. 5 I). IGV showed many hypomethylated DMRs in the gene body and CpGi in the promoter of FOXM1, which includes the CTCF binding site (Fig. S7C). Interestingly, this site has been reported to bind with SREBF2 (Fig. S7C). To assess whether FOXM1 regulates the core/marker genes of GICs, we used FDI-6 as a chemical inhibitor of FOXM1. FDI-6 down-regulated SOX4, SOX2, and PROM1 (Fig. 5J). This observation suggests that methionine metabolism altered the pluripotency of GICs through the cholesterol metabolism to the FOXM1 axis.

Methionine also targets protein synthesis
hsa-miR-33A has a DNA locus presenting many DMRs, adjacent to the SREBF2 locus (Fig. S7B), and is known to cause post-translational modifications of genes involved in cholesterol metabolism [7], so we analyzed non-coding RNAs using miRNA arrays. Exhaustive analysis of methionine depletion and statin administration identified one common SNORA17B (ACA43), but not hsa-miR-33A (Fig. 6A-D). ACA43 is predicted to induce pseudouridylation of U4938 of a 28S rRNA. Most pseudouridylation sites are located within functionally important ribosomal domains, which may affect the functional features of the ribosome. We thought that another common target of methionine and cholesterol metabolism might be a ribosomal protein. Entry of the mRNA expression data into GSEA using ontology gene sets found that gene ontology cytosolic large ribosomal subunit gene expression was extremely upregulated in an environment with adequate methionine (Fig. 6E). Identified RPL genes such as RPL26, RPL10A, RPL11, RPL9, and RPL34 are all constituent proteins of a large subunit of 28S, with the function to bind tRNA during the synthesis of polypeptides. We also confirmed that methionine depletion reduces protein synthesis with OPP assay (Fig. 6F).
Since epitranscriptome modifications such as pseudourinalization include methylation, we searched for genes involved in epitranscriptome methylation affected by methionine metabolic reorganization. The transcriptome Venn diagram showed methyltransferase-like protein 16 (METTL16), which is related to MAT2A, and methyltransferase-like protein 20 (METTL20) were predominantly up-regulated by methionine depletion (Table 1).

Discussion
SAM is the key molecule of transmethylation reactions and polyamine biosynthesis using methionine [8]. The biosynthesis of SAM is catalyzed by three major MATs, I, II, and III [9]. MATI and MATIII are encoded by MAT1A. MATII consists of an α2 catalytic subunit encoded by MAT2A and the β-regulatory subunit encoded by MAT2B. MATII functions as a transcriptional corepressor in the oxidative stress response, forming a SAM-integrating transcription regulation module that influences histone methyltransferase activity [10,11]. We demonstrated that methionine depletion caused decreased SAM and increased MAT2A and MAT2B in time-dependent manner. As a result, some genes are hypomethylated and others are hypermethylated, although the genome as a whole tends towards hypomethylation. Analysis of signaling pathways to determine what groups of genes are altered also shows a clear distinction between hyper-and hypo-methylation for each pathway. We speculate that this is probably where the reorganisation takes place, with the intention of gene expression to escape methionine depletion and allow the cell to survive. Although we did not search for the mechanism in this study, we believe that a search for the methylation status of histones and gene expression will be necessary in the future, as there are reports that SAM metabolism affects genome-wide histone methylation and is involved in the maintenance of pluripotent stem cells [12]. Since we expected decreased SAM to result in genomic hypomethylation of the promoter and gene body of MAT2A, we entered the RRBS data into IGV. Epigenetic methylation status of MAT2A and MAT2B were not so changed in spite of the marked increase in mRNA gene expression (Fig. S8). Thus, we searched for genes that affect the epi-transcriptome modulation using microarray results to understand the mechanism that causes the increase in key molecule MAT2A under the methionine-depleted state. We found two up-regulated N6-adenosine-methyltransferases, METTL16 and METTL20, which are involved in the post-transcriptional methylation of internal adenosine residues in eukaryotic mRNAs to form N6-methyladenosine). In agreement with our data, increased METTL16 occupancy on hp1 increased MAT2A splicing under SAM-limiting conditions [13]. Cells induce MAT2A expression by enhanced splicing of a retained intron. Interestingly, METTL16 also targets ACA43, one of the key factors shared by cholesterol and methionine metabolism, as discussed below [14].
Another up-regulated N6-methyladenosine methyltransferase, METTL20 (also known as electron transfer flavoprotein subunit beta lysine methyltransferase), is a mitochondrial lysine protein methylase that modifies Lys-199 and Lys-202 of the electron transfer flavoprotein from GSEA are shown on the right. Red column corresponds to upregulated genes and blue column corresponds to downregulated genes. NES: normalized enrichment score, FDR-q: false discovery ratio. F Methionine depletion reduced protein synthesis. We can detect the 5 FAM-Azide labeled protein in only cultured with methionine contained medium subunit beta, which is a mobile electron carrier coupling mitochondrial fatty acid oxidation to respiration [15,16]. Investigation using the knockout mouse showed METTL20 is a positive regulator of beta-oxidation. Betaoxidation is known as a metabolic node involved in glioma-genesis [17], and also provides an alternate source of ATP in nutrient unfavorable conditions. Our observation about METTL20 may represent the steps to obtain ATP.
Acetyl CoA is an important intermediate in the metabolism of lipids for energy storage and supply, from which fatty acids and cholesterol are synthesized. As mentioned above, beta-oxidation of fatty acids may be affected by methionine deprivation, but the metabolism of cholesterol was also greatly affected. Methionine depletion downregulates cholesterol biosynthesis and increases excretion. Furthermore, methionine depletion and statin administration showed similar phenotypic changes in GICs. These observations suggest that cholesterol metabolism is necessary to maintain the pluripotency of GICs. Biosynthesis of cholesterol is tightly regulated by feedback pathways because cholesterol is critical for cell growth and function. Cholesterol is important as a component of the plasma membrane and lipid rafts, and as a precursor for steroid hormones, bile acids, and vitamin D [18,19]. Other organs are more dependent on cholesterol supply from the bloodstream, whereas the brain obtains cholesterol exclusively from de novo synthesis, which involves the mevalonate and Bloch and Kandutsch-Russell pathways [20][21][22].
Within a tumor mass, imbalances in methionine concentration as well as oxygen concentration should exist. Using the Ivy GAP database, we compared the expression of genes involved in cholesterogenesis in the pseudopalisading cells and infiltrating tumor portions [23]. Our finding that cholesterogenesis was higher in the infiltrating tumor than in the pseudopalisading cell areas was consistent with the results of in vitro experiments. De novo synthesis is mediated by the transcriptional activity of sterol regulatory element binding proteins (SREBPs), which promote the transcription of enzymes involved in the cholesterol and fatty acid biosynthetic pathways, such as FDPS and HMGCR . When cholesterol levels fall below homeostatic levels, SREBP cleavage activating protein is isolated from insulin-inducible genes (INSIG1, INSIG2) resident in the endoplasmic reticulum [18,24]. The SREBP is then transported to the Golgi, where the active site is cleaved and activated as a transcription factor [25]. Among the three types (SREBP1a, SREBP1c, and SREBP2), SREBP2, which is encoded by the SREBF2 gene, mainly regulates genes involved with cholesterol biosynthesis. SREBF2 is upregulated with statin treatment by feedback but down-regulated by methionine depletion Table 1 Genes significantly up-regulated with change cut-off at 2.0 fold and P value cut-off at 0.05 in all MZGC1, MZGC2, and MZGC3 during methionine depletion METTL16 and METTL20 are depicted in red ( Fig. 5G and J). Such downregulation of SREBF2 under methionine depletion is thought to be the main cause of the downregulation of SREBF target genes of cholesterol biosynthesis.
Other groups have evaluated the effects of SREBPs on glioblastoma multiform development as a key transcription factor in the regulation of sterol homeostasis [26,27]. Under hypoxia and serum deprivation conditions, SREBP is upregulated to maintain the expression of fatty acid and cholesterol biosynthesis genes in glioblastoma multiform cells, and inhibition of SREBP activity under hypoxia results in the death of glioblastoma multiform cells [18,27]. So how does SREBF2 decrease under methionine depletion? We cannot determine the detailed mechanism, but the gene body of SREBF2, especially the 3' side including the hsa-miR-33A locus, is dramatically hypomethylated, and this might be the cause of downregulation of SREBF2. The hsa-miR-33A is also known as regulator of cholesterol efflux through ABCA1 gene expression [28].
To identify which miRNA contributes to the association between cholesterol metabolism and methionine depletion, we decided to explore the commonalities of changes in the miRNA array with methionine depletion and statin administration. Unexpectedly, only ACA43 but not hsa-miR-33A was downregulated in both arrays. ACA43 belongs to the H/ACA box class of small nucleolar RNAs based on the predicted hairpin-hinge-hairpintail structure and a conserved H/ACA-box motif. ACA43 is predicted to induce pseudouridylation of U4938, a 28S rRNA. Pseudouridylation is the isomerization of the nucleoside uridine to a different isomer, pseudouridine. Interestingly, ACA43 was originally cloned from HeLa cells [29] and is associated with the GAR1 protein, one of the components of the telomerase complex [30]. Taken together, these reports and our findings may indicate that cholesterol metabolism and methionine metabolism are also involved in cell immortalization via ACA43. At least, cholesterol and methionine metabolism altered the translation of mRNA through downregulated ACA43. Translation was also suppressed through down-regulation of the large subunit of ribosomal protein 28S by methionine depletion (Fig. 6E and F). The novel anticancer drug rocaglamide A, which inhibits translation by binding to the translation initiation factor eIF4A, inhibits the proliferation of target cancer cells, suggesting that translation inhibition due to methionine depletion may have inhibited the growth of GICs [31,32].
The final question is whether cholesterol metabolism is related to the apoptosis and pluripotency of GICs. Pharmacological inhibitors acting downstream of the mevalonate pathway induce cell death in glioma and pluripotent stem cells [33,34]. Inhibition of FDPS by pharmacological inhibitors and small interfering RNA alters the self-renewal of glioma stem cells [35]. HMGCS1 contributes to gastric cancer progression through activating Oct4 and SOX2 promoters [36]. These observations support our findings that methionine depletion promotes down-regulation of SOX2, OLIG2, PROM1, and FOXM1 along with down-regulation of cholesterol de novo biosynthesis. FOXM1 had a positive correlation with mevalonate pathway-related genes including HMGCR and SREBP2 in patients with hepatocellular carcinoma. Moreover, knockdown of HMGCR reduced FOXM1 expression in hepatocellular carcinoma [6]. We observed FOXM1 was increased by statin administration for a short period and decreased by methionine depletion along with SREBF2. Moreover, IGV showed first hypo-methylated gene body and CpGi of the FOXM1 promoter by methionine depletion and the SREBF2 binding site in CpGi of the FOXM1 promoter. These findings strongly suggest that methionine down-regulates FOXM1 through SREBF2. Since FDPS and HMGCS1 are targets of SREBF2, SREBF2 might be a key molecule-connecting cholesterol metabolism and pluripotency of stem cells. Finally, we found that the chemical inhibitor of FOXM1 downregulates core and marker genes of GICs. Taken together, we thought SREBF2-FOXM1 are important regulators of the core/marker genes in GICs (Fig. S9).
In conclusion, we found that methionine metabolism is closely related to cholesterol metabolism and ribosomal function, resulting in alteration of pluripotency, proliferation, and avoidance of apoptosis in GICs.
genes in the hyper/hypo-methylated categories by combining with the R package clusterProfiler. Supplementary Figure S4. Z-score of gene expression from Ivy GAP. Gene expression related de novo synthesis of cholesterol was upregulated in infiltrating tumor cells compared with pseudopalisading cells around necrosis. In contrast, the gene expression pattern related to efflux of cholesterol was exactly the opposite. Supplementary Figure S6. A qPCR of stem cell markers in MZGC4 and MZGC5. Statin treatment for 120 h down-regulated mRNA expression of FOXM1, PROM1, OLIG2,and SOX2. B Immunoblotting analysis showed that statin treatment for 2 days reduced SOX2protein level as well as methionine depletion for 7 days in MZGC7 and 8.