Nicotinamide riboside improves muscle mitochondrial biogenesis, satellite cell differentiation and gut microbiota composition in a twin study

Study participants

This study was approved by the Ethics Committee of the Helsinki University Central Hospital (protocol number 270/13/01/2008) and the study was conducted according to the principles of the Declaration of Helsinki. The study was registered at clinicaltrials.gov with the entry NCT03951285. Written informed consent was obtained from all subjects. Twenty-six MZ twin pairs belonging to the Finnish Twin Cohort study were screened for eligibility for the current study. Twenty-two MZ twin pairs with the cotwins having a BMI difference of at least 2.5 kg/m² were considered as BMI-discordant. Two such twin pairs did not meet the inclusion criteria, as shown in Figure S1. Another four twin pairs were considered as BMI-concordant (within-pair difference in BMI <2.5 kg/m²). The twins were 29-65 years old, with a mean age of 40. Fourteen of the MZ twin pairs were female and twelve were male pairs. The study participants had the following medications: hypertension (n=2; angiotensin II receptor blocker, calcium channel blocker, diuretic), asthma or allergy (n=7; an inhaled corticosteroid, beta-2-agonist, antihistamine), type 2 diabetes (n=1; metformin, insulin), rheumatoid arthritis (n=2; hydroxychloroquine, methotrexate), psychiatric (n=3; selective serotonin intake inhibitor, valproate, a tricyclic antidepressant, neuroleptic), oral contraceptives (n=2). The number of smokers is shown in Table S1.

Inclusion criteria were as follows: (i) age >18 years; (ii) BMI >18.5 kg/m² in both members of the twin pair; (iii) agreement to maintain current level of physical activity throughout the study; (iv) agreement to avoid vitamin supplementation or nutritional products with vitamin B3 14 days prior to the enrolment and during the study.

The exclusion criteria were: (i) unstable medical conditions as determined by the investigator; (ii) clinically significant abnormal lab results at screening (e.g. aspartate aminotransferase and/or ALT > 2 × upper limit of normal, and/or bilirubin > 2 × upper limit of normal); (iii) subjects who would have had a planned surgery during the course of the trial; (iv) History of or a current diagnosis of any cancer (except for successfully treated basal cell carcinoma diagnosed less than 5 years prior to screening). Subjects with cancer in full remission more than 5 years after diagnosis were acceptable; (v) history of blood/bleeding disorders; (vi) immunocompromised individuals such as subjects that had undergone organ transplantation or subjects diagnosed with human immunodeficiency virus; (vii) hepatitis; (viii) blood donation in the previous 2 months.

Subject recruitment and study design

BMI-discordant and BMI-concordant MZ twin pairs were screened from two large population-based cohorts of Finnish twins (FinnTwin12, Finntwin16, n=5000, 1200 of which MZ), established by Professor J.K. (Kaprio, 2013) and recruited by K.H.P. The procedures for the selection of study participants are described in Figure S1. One-to-one phone conversations were conducted at least once a month to follow up compliance to the study and potential symptoms or side effects. All participants were instructed to continue their normal routine and not make any changes to their physical activity and diet.

The examination protocol of this non-randomized, open-label study is described in Figure 1A. All twins from the BMI-discordant pairs received NR supplementation. For BMI-concordant twin pairs, one cotwin received NR and the other cotwin received placebo. The cotwins were allocated to the two groups randomly and in a blinded fashion. The NR dose, 1 g/day, was selected based on a previous study showing that 1 g/day niacin (another form of vitamin B3) supplementation had a favorable effect on systemic NAD⁺ levels, body composition and mitochondrial biogenesis (Pirinen et al., 2020). The weekly dose was initially 250 mg/week. The daily NR dose was gradually escalated by 250 mg/week so that the full dose of 1 g/day was reached in one month. The supplementation was continued for four months with the full dose. At the end of the study, the dose was slowly decreased to the initial level at a 250 mg/week rate. The collection of samples such as blood, WAT and muscle biopsies, saliva and feces was performed at the baseline and at the end of the study (i.e. 5 months). Four BMI-discordant twin pairs discontinued the intervention as described in Figure S1. Sixteen BMI-discordant and four BMI-concordant twin pairs were included in the further analyses, as shown in Figure S1. Safety laboratory tests included blood count, iron levels, and liver and kidney function tests. The analysis of body composition, clinical variables and whole blood NAD⁺ metabolites was performed for all study participants. The molecular measurements including gene expression, imaging and omics analyses were conducted only for the samples from the BMI-discordant twin pairs.

Body Composition, Energy Intake and Expenditure

Weight, height, waist-circumference, whole-body fat (assessed by using dual-energy x-ray absorptiometric scans), abdominal subcutaneous and visceral fat (by magnetic resonance imaging, MRI), and liver fat (by magnetic resonance spectroscopy, MRS) were measured as described previously in all study participants (Graner et al., 2012). In short, MRI/MRS for amounts of subcutaneous adipose tissue, visceral adipose tissue, and liver fat were performed on a clinical 1.5-T imager (Avanto, Siemens, Erlangen, Germany). To allow measurement of abdominal fat distribution, a stack of abdominal T1-weighted MR images (16 slices, slice thickness 10 mm, TR of 91 ms, TE of 5.2 ms and flip angle of 80°) were obtained from 8 cm above to 8 cm below the L4/5 lumbar intervertebral disks using frequency selective fat excitation. The amount of visceral and subcutaneous adipose tissue was determined from each slice using segmentation software SliceOmatic 5.0 (TomoVision, Quebec, Canada) utilizing region growing routine.

In turn, a point resolved spectroscopy sequence was used for volume selection in hepatic MRS. A 25×25×25 mm³ voxel was placed in the middle of the right liver lobe and liver spectra with TE of 30 ms and four averages were collected. Signal acquisition was triggered to end exhalation using a navigator belt to eliminate motion artifacts due to respiratory motion so that TR was kept >4000 ms. Liver spectra were analyzed with jMRUI 6.0 software (Stefan et al., 2009) and intensities of methylene and water resonances were determined using the AMARES algorithm (Vanhamme et al., 1997). Signal intensities were corrected for relaxation effects and liver fat was determined as an intensity ratio of methylene/(methylene+water). Ratios were further converted to mass fractions as described previously by (Kotronen et al., 2009).

Muscle mass and energy expenditure were measured via bioelectrical impedance analysis using the Tanita MC-980 device. Physical activity was measured using the Baecke questionnaire (Baecke et al., 1982) and dietary intake was assessed from 3-day food records and analyzed with the Diet32 program (Aivo Finland Oy), based on a national Finnish database for food composition (Fineli, www.fineli.fi).

Blood Laboratory Examinations

Routine blood tests were performed for all study participants. Blood samples were collected after overnight fasting, and whole blood, separated plasma and serum samples were frozen at −80C. Blood count, ALT, aspartate aminotransferase, HbA1c, glucose, insulin, C-peptide, lactate, total cholesterol, HDL, LDL, triglycerides, homocysteine and high-sensitivity CRP were analyzed using standardized methods at the HUSLAB laboratories.

A 75-g oral glucose tolerance test with 4-time points (0, 30, 60, and 120 min) was performed after a 12-hour overnight fast, followed by measurements of plasma glucose with spectrophotometric hexokinase and glucose-6-phosphate dehydrogenase assay (Roche Diagnostics) and serum insulin with time-resolved immunofluorometric assay (PerkinElmer). Plasma HDL was measured as cholesterol after precipitation of other lipoproteins with heparin–Mn2⁺ (Bachorik and Albers, 1986). Free fatty acids were determined using NEFA kit (Wako Chemicals), apolipoprotein B levels were determined with Apolipoprotein B Konelab Kit (Thermofisher Scientific), while plasma adiponectin levels were measured using ELISA kit (R&D Systems).

Adipose Tissue Sampling

Adipose tissue biopsies were collected from all study participants at baseline and at the end of the study (Figure 1A). The subcutaneous WAT biopsies were taken from superficial abdominal adipose tissue near the umbilicus using a needle biopsy under local anesthesia (lidocaine). WAT samples for RNA, DNA and metabolite analyses were snap frozen and stored in liquid nitrogen or −80°C until used.

Samples used for determination of the adipocyte size and number were, in turn, immediately processed by collagen digestion and separation of adipocytes by centrifugation as described previously by (Heinonen et al., 2014). In brief, the subcutaneous adipose tissue was minced and incubated for 1 h at 37 °C with continuous shaking in 10 ml of adipocyte medium (DMEM/F-12 (1:1) (Invitrogen) with 16 μmol/l biotin, 18 μmol/l panthotenate, 100 μmol/l ascorbate and antibiotic-antimycotic (Invitrogen)), supplemented with 2% bovine serum albumin (Sigma) and with 2 mg/ml collagenase A (Roche Diagnostics). Digestion was terminated by adding adipocyte medium supplemented with 10% newborn calf serum (Sigma) and centrifuging for 10 min at 600 × g. After washing the adipocytes with adipocyte medium, images of the adipocytes were taken with a light microscope (Zeiss, Axioplan2) using 50X magnification. Adipocyte diameters were automatically measured from the images using a custom algorithm for ImageJ (ImageJ 1.42q/ Java 1. 6.0 10 32-bit; https://github.com/birgittavdkolk/vanderkolk_etal_2021) (van der Kolk et al., 2021) which preprocessed the image to enhance the borders of the adipocytes and then used a circle-detection algorithm to identify the cells. The algorithm was tuned to identify the adipocytes taken using the standardized microscope settings and validated against 2000 manually measured diameters from 20 pictures (r = 0.85, p < 0.001). Mean adipocyte volume was calculated for each individual using the following formula: V = cell volume (μm³), d = cell diameter (μm), w=cell weight (μg), n = total adipocyte number, where m = total body-fat mass (kg) and V (m³)*rho (density, fat) = mean weight of a single adipocyte. Density of fat cell triglycerides = 0.915 g/ml. Adipocytes were assumed to be spheres.

Muscle Sampling, TEM, and Satellite Cell Histology and Quantification

Muscle tissue biopsies were collected at baseline and at the end of the study from all study participants (Figure 1A). A Bergström needle biopsy from the vastus lateralis muscle was taken under local anesthesia in sterile conditions by making an incision through the skin and underlying fascia. The muscle samples for RNA, DNA and metabolite analyses were snap frozen and stored in liquid nitrogen or at −80°C.

Samples for TEM analysis were fixed in 2.5% glutaraldehyde. For plastic embedding, they were then treated with 1% osmium tetroxide dehydrated in ethanol and embedded in epoxy resin. One mm section was stained with methylene blue (0.5% w/v) and boric acid (1% w/v). These methylene blue sections were used to mark the interesting areas for further TEM analyses by examining the sections with a light microscope. Thereafter, ultrathin (60– 90 nm) sections were cut on grids and stained with uranyl acetate and lead citrate by the Viikki Electron Microscopy Unit of the Institute of Biotechnology (EMBI, Biocenter Finland) and viewed with JEM-1400 transmission electron microscope (Jeol). From each study subject, 15 TEM images (both pre and post NR) were taken from the intermyofibrillar area of type I muscle fibers, determined by mitochondria content, using a 2000X magnification. Total mitochondria number per 10 µm² relative to the muscle fiber area, mitochondrial surface area per muscle fiber area (%), dimension values (average diameter (µm) and perimeter (µm)), and total muscle fiber area (µm²) were determined by painting each mitochondrion in Microscopy Image Browser (MIB) (Belevich et al., 2016). Further, these values were used to calculate the total mitochondria number per 10 µm² relative to the muscle fiber area, mitochondrial surface area per muscle fiber area (%, mitochondria density), and major shape descriptors, branching (form factor) and length-to-width ratio (aspect ratio).

Samples for muscle histological stainings were, in turn, snap frozen in liquid nitrogen-cooled isopentane. Cryosections from these samples were used to measure PAX7 immunostaining in muscle tissues: 9 µm thick frozen sections were post-fixed in 4% paraformaldehyde in DPBS for 10 min at RT and washed three times for 5 min in DPBS. Slides were then incubated in primary antibody (mouse monoclonal Pax7, Developmental Studies Hybridoma Bank (DSHB, Pax7 was deposited to DSHB by Kawakami A), 1:100 dilution in blocking solution (15% normal goat serum, 1.5% BSA, 0.5% Triton X-100 in 1X DPBS) overnight at 4°C. Thereafter, slides were washed three times for 10 min in DPBS and incubated in secondary antibody (Cy3 goat anti-mouse IgG1, 1:500 dilution, Jackson ImmunoResearch) and 1 μg/ml Hoechst for 1h at RT in order to visualize nuclei. Slides were washed as described above and mounted. Images were generated using Pannoramic 250 FLASH II digital slide scanner (3DHistec, Thermofisher Scientific) and Leica SP8 STED confocal microscope (Leica). For PAX7+ cell quantification, muscle sections were scored from a minimum of 2-3 sections/sample/condition and normalized to cross-sectional area of the muscle section scored. All quantifications were performed using Case Viewer 2.0 software (3DHistech, Thermofisher Scientific) and images were blinded for quantification. Details of the antibodies used are listed in the Key resources table.

Myoblast Cell Cultures

Muscle samples for primary muscle cell cultures, i.e. myoblast cultures were processed immediately after excision. In short, samples were washed with 5 ml isolation medium (0.05% trypsin-ethylenediamine tetraacetic acid (EDTA) (SAFC Biosciences), 1.8 g/l glucose (Sigma), 7.15 g/l HEPES (Sigma), 7.6 g/l NaCl (Sigma), 0.224 g/l KCl (Sigma), 1.2 mg/l phenol red (Sigma)). Thereafter, 5 ml isolation medium was added, and tissue was minced. Another 10 ml of isolation medium was added and the sample was transferred to a 20 ml beaker with a magnetic stirrer and placed in a 37°C incubator for 20 min. The supernatant containing muscle cells was transferred into a 50 ml centrifuge tube containing 15 ml of myoblast growth medium (Ham’s F10 with Glutamax (Gibco) supplemented with 15% fetal bovine serum (FBS, Thermofisher Scientific), 500 mg/l fetuin (Sigma), 200 mg/l insulin (Sigma), 0.05% bovine serum albumin (BSA, Sigma), 0.39 mg/l dexamethasone (DMSO) (Sigma), 50 mg/l gentamicin (Sigma), 2 mg/l human epidermal growth factor (Thermofisher Scientific)). This process was repeated a further two times to break down any remaining tissue. Finally, the supernatant containing the isolated muscle cells was centrifuged for 10 min at 100 × g to collect the cells and the resulting cell pellet was resuspended in 5 ml of myoblast growth medium. Myoblasts were first seeded on 25 cm² flasks (Corning). After reaching 70% confluence, they were subcultured in 150 cm² flasks (Corning) and finally on six-well plates (Corning) at a density of 100000 cells per well. When confluent, they were collected for RNA isolation and stored at −80°C.

Targeted quantitative NAD metabolome analysis

The measurement of whole blood NAD⁺ metabolome was carried out from all study participants after dual extractions as recently described (Pirinen et al., 2020). For analysis of Me4PY (group A analyte), samples were spiked with 400 pmol of [d3]-Me4PY (internal standard A). For analysis of NAD+, NADP, NMN, NAR, NAAD and ADPR, samples were dosed with 13C-yeast extract (internal standard B) as described (Elhassan et al., 2019; Trammell and Brenner, 2013).

In short, to quantify Me4PY (group A analyte), 75 μl of whole blood was added to 20 μL group A internal standard and 500 μl of 3:14%trichloroacetic acid (TCA):acetonitrile was added. The mixture was allowed to sit on ice for 20 min, after which samples were sonicated twice for 20 seconds and centrifuged at 4°C for 13 minutes at 16 100 × g. Next, the supernatant was removed and dried under vacuum overnight at room temperature. The samples were reconstituted in 2% acetonitrile/water immediately prior to the analytical run. To quantify NAD⁺, NADP, NMN, NAR, NAAD and ADPR (group B analytes), 75 μl of whole blood was added to 20 μl group B internal standard prepared in water and mixed with 500 μl of 3:1 4% (TCA):acetonitrile with vortexing. After resting on ice, the samples were centrifuged as described above. Next supernatant was removed and dried under vacuum overnight, and reconstituted in 2% acetonitrile/water. After reconstitution, samples were transferred to Waters polypropylene plastic total recovery vials and stored in a Waters Acquity H class autosampler maintained at 8°C until injection. In all cases, 8 μl of extract was loaded onto the column. For one set of samples, group B analytes plus Me4PY were extracted as for the group B analytes. 100 μl of blood was used and 200 pmol d3-Me4PY was added to the internal standard mix. All analytes in this group except NAAD had a corresponding stable labeled internal standard. Labeled NAD⁺ was used as the internal standard for NAAD.

Separation and quantitation of analytes were performed with a Waters Acquity LC interfaced with a Waters TQD mass spectrometer operated in positive ion multiple reaction monitoring mode as described (Trammell and Brenner, 2013) with minor modifications. MRM transitions monitored were: NAD⁺ 664.1>136.1; 13C10-NAD⁺ 674.1>136.1; NADP 744.1>136.1; 13C10-NADP 754.1>136.1; NAAD 665.1>136.1; NMN 335>123: 13C5-NMN 340>123; NAR 256>124; 13C5-NAR 261>124; ADPR 560>136.1; 13C10-ADPR 570>136.1; Me-4-py 153>136; d3-Me-4-py 156>139. Separate Hypercarb columns (Thermo Scientific) were used for the separations. The conditions for group A were solvent A 10 mM NH4OAc with 0.1% formic acid, solvent B acetonitrile with 0.1% formic acid, solvent D methanol, flow 0.30 mL/min; gradient initial 98% A, 2% B; 2.25 min, 98% A 2%B; 11 min, 74.8% A 17.9% B 7.3% D; 11.1 min,10% A 90% B; 14.3 min 10% A 90% B; 14.4 min, 98.2% A 2%B; end 18.5 min. Conditions for Group B were solvent A 7.5mM NH4OAc with 0.05% NH4OH, solvent B acetonitrile with 0.05% NH4OH; flow 0.353 ml/min; gradient initial 97% A 3% B; 1.8 min, 97% A 3% B; 10 min 65.5% A 34.5% B; 11 min 10% A 90% B; 13.2 min, 10% A 90% B; 13.3 min, 97% A 3% B; end 19 min.

mtDNA Analyses and Quantification

Total DNA, including mtDNA, was extracted from the snap-frozen muscle and WAT samples (15 mg and 150 mg, respectively) using standard phenol–chloroform extraction and ethanol precipitation method, and stored at +4°C until use. mtDNA content was determined by analyzing the ratio of mtDNA amount relative to genomic DNA using SYBRGreen-based quantitative PCR (qPCR). To detect genomic DNA, primers for nuclear amyloid β precursor protein (APP), beta-2-microglobulin (B2M) and hemoglobin subunit beta (HBB) were used. In turn, to measure mtDNA, mitochondrial 16S rRNA (16S), cytochrome B (CYTB) and 12S rRNA (DLOOP) primers were chosen. The PCR primers are provided in the Key resources table.

qPCR was performed using 1X SYBRGreen MasterMix (Thermofisher Scientific), 10 µM of primers and 2 ng of DNA template in 384 well plate format using CFX384 Touch™ Real-Time PCR Detection Systems (Bio-Rad Laboratories). Thermal cycling as triplicates was run using initial denaturation of 3 min at 95°C, 39 cycles of 10 s at 95°C and 30 s at 62°C, final extension of 10 s at 95°C and melting curve analysis step from 65°C to 95°C (raised gradually in 0.5°C increments). The mtDNA amount was presented as a relative level of mtDNA genome per nuclear genome. Data analysis was performed using the ΔΔCt method with the qBASE+ 3.2 software (Biogazelle).’

Gene expression

Total RNA from myoblasts was isolated using TRIzol reagent (Invitrogen) according to the manufacturer’s instructions. In turn, total RNA from muscle tissue and WAT was isolated using AllPrep DNA/RNA/miRNA Universal Kit (Qiagen) according to the manufacturer’s instructions with a few minor modifications. In short, muscle tissue (30 mg in 600 µl Buffer RLT Plus) was lysed using Tissuelyser II (Qiagen) at 25 Hz for 2 min and the lysate was transferred into the AllPrep DNA Mini spin column. WAT (250 mg in 1 ml Buffer RLT Plus), in turn, was lysed using Tissuelyser II (Qiagen) at 25 Hz for 2 × 2 min. As the total sample volume was greater than recommended, the lysate was divided into two AllPrep DNA Mini spin columns. From this point on, muscle and WAT samples were processed identically to the manufacturer’s instructions. RNA was stored at −80 °C until use.

Complementary DNA (cDNA) was synthesized from 0.5 µg of muscle tissue and myoblast RNA using QuantiTect Reverse Transcription Kit (Qiagen) according to the manufacturer’s instructions. In turn, SuperScript VILO cDNA Synthesis Kit (Thermofisher Scientific) was used to synthesize cDNA from 0.5 µg of WAT RNA, according to the manufacturer’s instructions.

Expression levels of genes associated with NAD⁺ biosynthesis, mitochondrial metabolism, mitochondrial OXPHOS, and additional potential NR target genes were determined from tissue samples with SYBRGreen-based quantitative reverse transcription PCR (qRT-PCR). Additionally, gene expression levels of myogenic marker MYOG and stem cell marker PAX7 were measured from myoblasts in order to assess their differentiation potential. The PCR primers for these genes are provided in the Key resources table. qRT-PCR was performed using 1X SYBRGreen MasterMix (Thermofisher Scientific), 10 µM of primers and 5 ng of DNA template in 384 well plate format using CFX384 Touch™ Real-Time PCR Detection Systems (Bio-Rad Laboratories). Samples were run in duplicates with conditions otherwise identical to mtDNA quantification PCR. To normalize the expression data, three reference genes were used: actin, beta-2-microglobulin and tyrosine 3-monooxygenase/tryptophan 5-monooxygenase activation protein zeta. Data analysis was performed with the ΔΔCt method using the qBASE+ 3.2 software (Biogazelle).

DNA methylation analysis

High molecular weight DNA was extracted from WAT and muscle using QiAmp DNA Mini kit (Qiagen) and bisulfite converted by EZ DNA Methylation kit (ZYMO Research) following the manufacturer’s protocol. DNA methylation was quantified using Infinium MethylationEPIC BeadChip arrays (Illumina) according to the manufacturer’s instructions. The data was preprocessed using R-package minfi (Aryee et al., 2014). Sample quality control was performed in R-package MethylAid with default settings (van Iterson et al., 2014). Furthermore, bad quality probes were removed with the following criteria: (i) detection p-value > 0.01, (ii) bead counts < 3, (iii) zero methylation probe intensity and (iv) ambiguously mapping probes identified by Zhou and colleagues (Zhou et al., 2017). The data were normalized using functional normalization (Fortin et al., 2014) followed by Beta Mixture Quantile normalization (Teschendorff et al., 2013). M values representing the methylation status of each probe were calculated by taking the log2 ratio of methylated probe intensity to unmethylated probe intensity. CpG sites that annotated to genes from which there was gene expression data available in either of the tissues, were selected for the analysis, resulting in a final dataset of 518 probes in WAT and 619 in muscle. Names of the genes are listed in the Key resources table.

Methylation levels in DNA repetitive elements (e.g. LINE-1) across the genome can be used as a surrogate measure for the global DNA methylation level of a sample. These elements are often hypermethylated, in contrast to CpGs located in other genomic areas (Rollins et al., 2006). The average LINE-1 methylation levels were computed using R-package REMP (Zheng et al., 2017) for each individual. Paired Wilcoxon signed-rank test was applied to compare the global methylation levels pre versus post NR across all twins from the BMI-discordant pairs (regardless of heavy/lean status), as well to compare global methylation levels as a result of NR supplementation, i.e. delta variables (change from the baseline to 5 months). The level of significance (two-tailed) was set at p<0.05.

Differential methylation analysis was performed for the selected CpGs in WAT (n=14 twin pairs / 28 individuals) and muscle (n=12 twin pairs / 25 individuals) using a linear model (package limma in R-Bioconductor) to identify those CpG sites that showed significantly different methylation levels (Benjamini-Hochberg FDR<0.05) pre versus post NR intervention across all twins from the BMI-discordant pairs (regardless of heavy/lean status). Regression analysis was also performed to identify CpGs that had different methylation levels as a result of NR intervention, i.e. delta variables (change from the baseline to 5 months). The regression models were adjusted for leaner/heavier status, sex, age, smoking status and twinship.

Plasma metabolite isolation and non-targeted metabolite profiling

We performed, at the LC-MS metabolomics center (Biocenter Kuopio, University of Eastern Finland), non-targeted metabolite profiling on plasma metabolites. We analyzed samples using two different chromatographic-mass spectrometric techniques, which were reversed-phase (RP) and hydrophilic interaction (HILIC) chromatography techniques coupled online to high-resolution mass spectrometry, and further acquired data with positive and negative electrospray ionization (ESI) mode (Noerman et al., 2021; Pekkinen et al., 2013).

Briefly, plasma samples were thawed on ice and a 100 μl aliquot of plasma was dispensed into a 96-well Captiva ND filter plate 0.2 µm PP (Captiva ND, 0.2 µm PP, Agilent Technologies) containing 400 µl of ice-cold acetonitrile. Samples were mixed with a pipette to thoroughly precipitate plasma proteins. Samples were then centrifuged at 700 × g for 5 min at 4°C and the supernatants were collected to a 96-well storage plate and stored at −10°C.

A pooled sample from all biological samples per experiment was injected at the beginning of the sequence to equilibrate the analytical platform and then after every 12 samples throughout the analysis for quality control (QC). In addition, a solvent blank was prepared and injected at the beginning of the sequence. Samples were randomized before the analysis.

The analysis of amphiphilic metabolites was carried out using an ultra-high performance liquid chromatography (Vanquish Flex UHPLC system, Thermo Scientific, Bremen, Germany) coupled online to high-resolution mass spectrometry (HRMS, Q Exactive Classic, Thermo Scientific, Bremen, Germany). Samples were analyzed using the RP technique. The sample solution (2 µl) was injected onto an RP column (Zorbax Eclipse XDBC18, 2.1 × 100mm, 1.8 μm, Agilent Technologies, Palo Alto, CA, USA) that was kept at 40°C. The mobile phase, delivered at 400 μl/min, consisted of water (eluent A) and methanol (eluent B), both containing 0.1 % (v/v) of formic acid. The following gradient profile was used: 0–10 min: 2 to 100% B, 10–14.50 min: 100% B, 14.50–14.51 min: 100 to 2% B; 14.51–20 min: 2% B. The sample tray was at 10°C during these analyses. Mass spectrometry was equipped with heated electrospray ionization (ESI). The positive and negative ionization modes were used to acquire the data in centroid mode. At the beginning and end of the sample analysis, data-dependent product ion scans (MS2) were acquired for each mode. The detector was calibrated before the sample sequence and subsequently operated at high mass accuracy (<2 ppm). Continuous mass axis calibration was performed by monitoring reference ions m/z 214.08963 in the positive ionization mode.

For the analysis of hydrophilic compounds, an ultra-high performance liquid chromatography (1290 LC system, Agilent Technologies, Waldbronn, Karlsruhe, Germany) coupled online to high-resolution mass spectrometry (HRMS, 6540 UHD accurate-mass quadrupole-time-of-flight (qTOF) mass spectrometry, Agilent Technologies, Waldbronn, Karlsruhe, Germany) was used (Pekkinen et al., 2013). The sample solution (2 µl) was injected onto a column (HILIC, Acquity UPLC® BEH Amide 1.7 µm, 2.1×100 mm, Waters Corporation, Milford, MA, USA) that was kept at 45°C. Mobile phases, delivered at 600 μl/min, consisted of 50% (v/v)(eluent A) and 90% (v/v)(eluent B) acetonitrile, respectively, both containing 20 mM ammonium formate (pH 3). Following gradient profile was used: 0– 2.5 min: 100% B, 2.5–10 min: 100% B → 0% B, 10–10.1 min: 0% B → 100% B; 10.1–12.5 min: 100% B. The total runtime was 12.5 min. The sample tray was at 10°C during these analyses. Mass spectrometry was equipped with a heated ESI source, operated in both positive and negative ionization modes. For data acquisition, a 2 GHz extended dynamic range mode was used in both positive and negative ion modes from 50 to 1600 (m/z). Data was collected in the centroid mode. At the beginning and end of the sample analysis, data-dependent product ion scans (MS2) were acquired for each mode. The detector was calibrated before the sample sequence and subsequently operated at high accuracy (<2 ppm). Continuous mass axis calibration was performed by monitoring two reference ions from an infusion solution throughout the runs. The reference ions were m/z 121.050873 and m/z 922.009798 in the positive mode and m/z 112.985587 and m/z 966.000725 in the negative mode.

Further information about the instruments set-up and data acquisition parameters, in concert with molecular feature finding and peek picking can be obtained from the previous publications (Noerman et al., 2021; Pekkinen et al., 2013).

Differential global metabolite profile analysis and metabolite identification

Differential metabolite analysis was performed on the global metabolite levels of the twins from BMI-discordant pairs (n=14 twin pairs/28 individuals) using a linear model (package limma in R-Bioconductor) to identify metabolite features that had significantly different levels (nominal p<0.05) pre vs post NR across all study subjects (regardless of heavy/lean status). Regression models were adjusted for leaner/heavier status, sex, age, smoking status and twinship (shared genetic effects). Regression analysis was also carried out to identify metabolite features showing different levels in heavier compared to leaner individuals as a result of NR intervention, i.e. delta metabolites (change from the baseline to 5 months). The regression models were adjusted for sex, age, smoking status and twinship. Significant metabolites were then further annotated, and enriched pathways were identified using MetaboAnalyst (metaboanalyst.ca) and Ingenuity Pathway Analysis tool (IPA, Qiagen).

The final data set contained 580 and 2398 molecular features from HILIC and RP analysis, respectively. Annotation of these significantly different metabolites found between the study groups was generated based on accurate mass and isotope information; i.e., ratios, abundances, and spacing, as well as product ion spectra (MS2) against existing libraries, either in-house for level I or online spectral databases (i.e., Human Metabolome Database, MzCloud, METLIN, ChemSpider) for level II-IV, according to the guidelines from the Metabolomics Standard Initiative (Sumner et al., 2007).

Gut microbiota sequencing and composition analyses

Fecal samples were taken at home or at the hospital. For stool collection at home, study subjects were provided with the supplies and instructions for the collection. Study subjects were asked to keep the sample tube in a cold environment, and to bring the sample to the hospital. The stool was then immediately stored at −80 °C without any treatment. The interval between the passage of stool and storage at −80 °C was within 12 to 24 hours. The sample collection at the hospital was conducted in the same fashion. For the analyses of gut microbiota, the exclusion criteria were antibiotic course one month prior to the sample collection.

For the microbial community analysis, the total microbial DNA was extracted from ~80 mg of feces by first homogenizing the sample with a bead-beating method and then using Stool Extraction Kit v2 (Hain Lifescience GmbH, Nehren, Germany) and semi-automated GenoXtract (Hain Lifescience Gmbh). Thereafter, the rRNA gene was amplified using universal primers S-D-Bact-0341-b-S-17 and S-D-Bact-0785-a-A-21 targeting the V3-V4 regions of the SSU (small subunit) rRNA gene, details provided in the Key resource table. The PCR was performed using 1X SYBRGreen MasterMix (Thermofisher Scientific, Waltham, MA, USA), 0.5 µM of primers and 20 ng of DNA template. Thermal cycling for the first PCR was run using initial denaturation of 10 min at 95°C, 30 cycles of 30 s at 94°C and 60 s at 52°C, 60 s at 72°C and final extension for 5 min at 72°C. To add Ion Torrent PGM sequencing adapters and barcodes to the ends of the PCR product, one µl of the PCR product was used as a template in the second PCR, where 10 cycles were performed using linker and fusion primers (0.05 µM of M13_ S-D-Bact-0341-b-S-17, 0.5 µM of IonA_IonXpressBarcode_M13 and 0.5 µM P1_ S-D-Bact-0785-a-A-21), with conditions otherwise identical to the first amplification. Sequencing was performed using Ion Torrent PGM (Thermofisher Scientific). Prior to sequencing, the PCR products were purified with AMPure XP (Beckman Coulter), quantified with PicoGreen (Quant-iT PicoGreen dsDNA Assay Kit, Thermofisher Scientific) according to manufacturer’s instructions, and pooled in equimolar quantities for sequencing on Ion Torrent PGM using Thermofisher Scientific’s Hi-Q View OT2 Kit for emulsion PCR, Hi-Q View Sequencing Kit for the sequencing reaction, and Ion 318 Chip v2. The 16S rRNA gene sequences were quality-filtered and clustered to operational taxonomic units (OTUs) at 97% similarity using the CLC Microbial Genomics Package (Qiagen).

The rRNA gene sequences were classified using the SILVA SSU Ref database (99%, arb-silva.de/projects/ssu-ref-nr/). Statistical analyses for gut microbiota were performed with the CLC Microbial Genomics Package and SPSS Statistics (IBM). The alpha-diversity of the gut microbiota was quantified with the Shannon index and beta-diversity analysis was based on Bray-Curtis distance and PERMANOVA for significance testing. The taxonomic differences between pre versus post NR supplementation samples were analyzed with paired Wilcoxon signed-rank test in SPSS, followed by Benjamini-Hochberg FDR correction for multiple testing. The level of significance (two-tailed) was set at p<0.05 after the multiple testing correction.

Microbiota data were pre-processed and normalized using packages phyloseq and edgeR in R-Bioconductor. Differential analysis was performed for the gut microbiota levels of the twins (n=11 twin pairs/22 individuals) using a linear model (package limma in R-Bioconductor (Ritchie et al., 2015) to identify microbiota that had significantly different levels (nominal p<0.05) pre vs post NR across all study subjects (regardless of heavy/lean status). The regression models were adjusted for leaner/heavier status, sex, age, smoking status and twinship. Regression analysis was also carried out to identify differentially abundant microbiota showing different levels in heavier compared to leaner individuals as a result of NR intervention, i.e. delta microbiota (change from the baseline to 5 months). The regression models were adjusted for sex, age, smoking status and twinship.

F. prausnitzii qPCR analysis was performed using iQ SYBR Supermix and CFX96™ Real-time PCR Detection System (Bio-Rad Laboratories, Hercules, CA, USA). To amplify F. prausnitzii, 20 ng of template DNA and 10 pM of primers were used whereas to amplify the 16S rRNA gene, 10 ng of template DNA and 2.5 µM of primers were added per reaction. The PCR primers are provided in the Key resources table. The qPCR program for F. prausnitzii was as follows: pre-incubation at 95°C for 10 min, 40 cycles at 95°C for 30 sec, 60°C for 30 sec and 72°C for 15 sec. The qPCR program for the 16S rRNA genomic region was pre-incubation at 95°C for 3 min, 40 cycles at 95°C for 30 sec, 53°C for 30 sec and 72°C for 30 sec. The results of F. prausnitzii were normalized to the results of 16S rRNA, and the fold change between the pre and post NR samples was calculated with the ΔΔCt method. Due to a great interindividual variation, the effect of NR on F. prausnitzii amount was analyzed with Cohen’s d as follows: where M₁ and M₂ are the means for the 1st and 2nd samples, and SD_pooled is the pooled standard deviation for the samples.

Association analyses

Associations between CpG methylation (outcome variable) and gene expression (explanatory variables) were identified in study participants with the respective data available. First, delta gene expression (the change from baseline to 5 months) was calculated which was then compared to corresponding delta DNA methylation (the change from baseline to 5 months) of the given gene. Linear modelling was used to test for the associations, adjusting for age, sex, smoking status and twinship. Benjamini-Hochberg FDR correction was performed separately for each CpG-variable comparison with a significance threshold of <0.05.

Furthermore, the associations between NAD/global metabolites (outcome variable) and microbiota (explanatory variables) were examined. First, regression analysis was performed to identify changes, between 0 and 5 months, in NAD metabolites in relation to changes in the microbiota. Next, we identified, between 0 and 5 months, changes in NAD metabolites in relation to the baseline microbiota. Finally, a similar analysis was performed again, to identify changes, between 0 and 5 months, in global metabolites as well as clinical parameters, in relation to changes in the microbiota. Since our focus was on the metabolites that had already been determined to be differentially abundant in pre vs post NR samples, only their relationships with microbiota were examined. For relationships, i.e. associations, that were significant (FDR<0.05), pathway analysis was performed using the IPA tool. All of the regression models mentioned above were adjusted for leaner/heavier status, sex, age, antibiotics usage, smoking status and twinship.