Materials and methods

Dataset and experimental design

The dataset analyzed here was GSE85718, generated from the long-term NMN study by Mills et al. [6] and deposited in the Gene Expression Omnibus (GEO) [8]. The analysis was a secondary reanalysis of processed public data and did not involve new animal experiments, treatment assignments, sample collection, or wet-laboratory measurements. The dataset included 48 samples across three tissues: skeletal muscle, liver, and white adipose tissue (WAT) (Figure 1). Within each tissue, samples represented two ages, 6 months and 12 months; two treatments, control and NMN; and four biological replicates per age-by-treatment cell. The microarray platform was GPL6885, the Illumina MouseRef-8 v2.0 expression beadchip.

Figure 1: Schematic of the analysis pipeline(a) The GSE85718 design comprises 48 microarray samples spanning three metabolic tissues, two ages, two treatments, and four biological replicates per cell. (b) Probe-level data are collapsed to gene level by retaining the highest-mean probe per symbol and then quantile-normalized to yield the matrix used for modeling. (c) Within each tissue, expression is fit to an age × treatment linear model; the age-by-treatment interaction is the central test, and rescue candidates require an age effect plus an opposite-signed interaction. (d) Candidate lists are sequentially filtered by interaction effect size, statistical confidence, cross-tissue robustness with direction consistency, and pathway context via MitoCarta overlap and preranked gene set enrichment analysis (GSEA) leading-edge analysis. The schematic depicts the computational workflow only; no results are shown, and all proposed mechanistic routes remain hypotheses requiring independent validation.Image created by the author using PowerPoint (Microsoft® Corp., Redmond, WA).

Expression processing and normalization

Probe-level expression data were downloaded from GEO using GEOparse. Probes were mapped to gene symbols using the GPL6885 platform annotation and collapsed to gene-level values by retaining the highest-mean probe for each symbol. This reduced 25,697 probes to 18,120 genes.

The analysis used the processed GEO VALUE field rather than raw intensity files. The values appeared to be on a compressed or log-like scale because the maximum observed value did not exceed the threshold used in the analysis script to identify linear-scale intensities; accordingly, no additional log2 transformation was applied. Per-sample quantile normalization was then performed. This preprocessing decision was based on the available processed matrix and was not supported by a separate array-density, MA-plot, or raw-intensity diagnostic analysis. The resulting normalized matrix was used for all tissue-specific modeling.

No arrays were excluded as outliers, and no array weighting, surrogate-variable adjustment, batch correction, or cell-composition adjustment was performed. These omissions should be considered when interpreting the small-sample interaction results.

Statistical model

For each tissue, gene expression was modeled as expression ~ age × treatment. Age was encoded as 0 for 6 months and 1 for 12 months. Treatment was encoded as 0 for control and 1 for NMN. The age coefficient therefore estimated the age effect within the control group. The treatment coefficient estimated the NMN effect at six months. The interaction coefficient estimated whether NMN changed the aging slope. This interaction term was the main coefficient of interest because it directly tests whether NMN modifies age-related transcriptional change.

Ordinary least squares (OLS) was used because the primary scientific question was the age-by-treatment interaction coefficient itself and because the small per-cell sample size (n = 4) already limits power. The primary analysis therefore did not apply empirical-Bayes variance moderation. This choice provides a straightforward coefficient-level interpretation but may result in less stable variance estimates than moderated microarray models such as limma. A limma-based sensitivity analysis was not performed in the present computational run; therefore, equivalence between the OLS results and empirical-Bayes-moderated results cannot be claimed. The absence of this sensitivity analysis is an additional limitation and supports treating the candidate lists as hypothesis-generating.

Definition of NMN-rescue candidates

An NMN-rescue candidate was defined by three criteria. First, the gene had to show an age effect in control mice. Second, the gene had to show an age-by-treatment interaction. Third, the interaction coefficient had to have the opposite sign from the control aging coefficient. This sign-reversal criterion distinguishes a general NMN effect from a pattern consistent with NMN-associated attenuation of age-related expression drift. Because no genes passed strict genome-wide false discovery rate (FDR) thresholds, the rescue lists were generated using relaxed nominal criteria: P_age, the P value for the age effect in control mice, <0.05; P_interaction, the P value for the age-by-treatment interaction, <0.05; and opposite sign between the age effect and interaction term. These uncorrected nominal thresholds create a substantial false-positive risk, particularly when many thousands of genes are tested. The relaxed lists were used for candidate generation only and not for confirmatory inference.

Multiple testing was addressed using the Benjamini-Hochberg FDR procedure, computed by tissue and statistical term [9]. The use of a linear model follows the broad framework of microarray differential-expression modeling, although the present pipeline used OLS rather than empirical-Bayes moderation [10,11].

Downstream prioritization

Downstream prioritization proceeded in several steps. First, the relaxed rescue lists were filtered by interaction effect size using an absolute interaction coefficient threshold of 0.8. Second, high-confidence candidates were defined as genes with P_interaction <0.01 and an absolute interaction coefficient of at least 0.8. These high-confidence labels indicate stronger candidates within the nominal-P framework; they do not indicate genome-wide statistical significance. Third, cross-tissue robustness was assessed by identifying genes rescued in at least two tissues. Fourth, direction consistency was evaluated by asking whether the interaction coefficients shared the same sign across tissues (Appendix 1).

Fifth, mitochondrial overlap was assessed using a MitoCarta-informed approach. Because the full mouse MitoCarta3.0 list was not available in the run environment, a curated fallback mitochondrial set of 66 symbols was used. The exact fallback list is embedded in the analysis script used for this study. Because this list is incomplete relative to the full MitoCarta3.0 mouse annotation, all mitochondrial overlap and enrichment findings are provisional. MitoCarta3.0 remains the preferred reference because it provides an updated mitochondrial proteome with sub-organelle localization and pathway annotations [12].

Gene set enrichment analysis (GSEA)

GSEA was performed using preranked interaction coefficients [13]. Genes were ranked by the raw, unscaled NMN-by-age interaction coefficient. No coefficient standardization or additional scaling was applied. Preranked GSEA used 1,000 permutations, a minimum gene-set size of 10, and a maximum gene-set size of 500. The ranking asked which pathways were associated with positive or negative NMN-by-age interaction effects.

Gene Ontology (GO) Biological Process 2021 and the Kyoto Encyclopedia of Genes and Genomes (KEGG) 2019 Mouse were treated as the more reliable mouse-oriented libraries. Hallmark gene sets were interpreted more cautiously because they are human-derived and required upper-casing mouse gene symbols as an approximate orthology proxy. Leading-edge overlap was then used to identify robust genes that appeared repeatedly within pathway-driving subsets. Leading-edge membership was used as a prioritization heuristic and not as independent evidence that the corresponding pathway or cellular process changed functionally.

Software and code availability

The pipeline was implemented in Python and used GEOparse for GEO retrieval, pandas and NumPy for data handling, SciPy for statistical calculations, statsmodels for Benjamini-Hochberg correction, scikit-learn for quantile normalization, matplotlib and seaborn for visualization, and GSEApy for preranked GSEA (Appendix 2). The primary model was implemented with vectorized ordinary least-squares calculations, and false-discovery correction used the statsmodels multiple-testing functions. Exact package-version records were not captured during the original computational run. The complete annotated analysis script and generated processing outputs should therefore be retained with the study files to support exact rerunning of the analysis.