Gene‐Metabolite Network Linked to Inhibited Bioenergetics in Association With Spaceflight‐Induced Loss of Male Mouse Quadriceps Muscle

ABSTRACT Prolonged residence of mice in spaceflight is a scientifically robust and ethically ratified model of muscle atrophy caused by continued unloading. Under the Rodent Research Program of the National Aeronautics and Space Administration (NASA), we assayed the large‐scale mRNA and metabolomic perturbations in the quadriceps of C57BL/6j male mice that lived in spaceflight (FLT) or on the ground (control or CTR) for approximately 4 weeks. The wet weights of the quadriceps were significantly reduced in FLT mice. Next‐generation sequencing and untargeted mass spectroscopic assays interrogated the gene‐metabolite landscape of the quadriceps. A majority of top‐ranked differentially suppressed genes in FLT encoded proteins from the myosin or troponin families, suggesting sarcomere alterations in space. Significantly enriched gene‐metabolite networks were found linked to sarcomeric integrity, immune fitness, and oxidative stress response; all inhibited in space as per in silico prediction. A significant loss of mitochondrial DNA copy numbers in FLT mice underlined the energy deprivation associated with spaceflight‐induced stress. This hypothesis was reinforced by the transcriptomic sequencing–metabolomics integrative analysis that showed inhibited networks related to protein, lipid, and carbohydrate metabolism, and adenosine triphosphate (ATP) synthesis and hydrolysis. Finally, we discovered important upstream regulators, which could be targeted for next‐generation therapeutic intervention for chronic disuse of the musculoskeletal system. © 2020 The Authors. Journal of Bone and Mineral Research published by American Society for Bone and Mineral Research.

At the time of biomolecular extraction, the quadriceps tissues were first weighed and then cryogenically homogenized using a Cryomill (Retsch GmbH, Germany). An aliquot of the homogenates was allotted for extraction RNA/DNA using the standard TRIzol (Invitrogen, Thermo Fisher Scientific, Wilmington, MA) method.
DNA extracted from the interphase of TRIzol solution was assayed to measure the copy number variation (CNV) of mitochondrial DNA (mtDNA) using the NovaQUANT™ Mouse Mitochondrial to Nuclear DNA Ratio Kit (Novagen, WI). Vendor provided PCR primer pairs were used to amplify two mitochondrial genes (trLEV and 12s RNA) in Applied Biosystems 7900HT real-time qPCR. CNV was calculated following the equation: The sequencing library preparation was carried out using 500 ng of total RNA input with the TruSeq Stranded mRNA Sample Preparation Kit (Illumina Inc, San Diego, CA) following the manufacturer's procedures. In brief, total RNA was fragmented using divalent cations under elevated temperature. Cleaved RNA fragments were then converted to first strand cDNA using reverse transcriptase SuperScript II (Thermo Fisher, Waltham, MA). Strand selection was done by replacing dTTP with dUTP followed by second strand cDNA synthesis using DNA Polymerase I and RNase H (Illumina Inc, San Diego, CA). A tail was added to the 3' ends of the blunt fragments. Unique indexing adapters were ligated to the ends of the ds cDNA fragments and the products were enriched with PCR and purified to create the final cDNA library. The quality of the libraries were determined using the on D1000 ScreenTape on the Agilent TapeStation 2200 (Agilent Technologies, Santa Clara, CA, USA). Libraries were quantified with real-time quantitative PCR using KAPA SYBR® FAST Universal qPCR Kit (Kapa Biosystems). Next, the libraries were normalized, pooled, and quantified using the dsDNA HS Kit on the Qubit 3.0 Fluorometer ; followed by sequencing 2x150 cycles paired-end reads on an Illumina HiSeq 4000 platform (Illumina, Inc., San Diego, CA) using the Illumina HiSeq 3000/4000 SBS kit following the manufacturer's instructions.
mRNA raw reads were quality filtered with Q score above 30 and adaptors were trimmed using Trimmomatic v0.38 and reads were further aligned using STAR aligner v2.6.0 against UCSC genome assembly for Mus Musculus mm10. Gene counts were further filtered using counts per million (CPM) cutoff greater than 0 across eight samples. Counts were normalized using the TMM normalization method using edgeR package v3.24.3. The library depth and count distribution was plotted in Figure S2. Differential expression was calculated by fitting a negative binomial generalized log-linear model to the read counts for each gene using the limma package v3.38.3.
The second aliquot of the homogenized quadriceps tissues was re-suspended in ice cold Phosphate-buffered saline (PBS) with a protease inhibitor tablet. The resuspension was centrifuged at 5000x g for 5 minutes and the supernatant was collected. The load of myosin protein was estimated by using the Abbexa (Cambridge, UK) Mouse Myosin ELISA kit (Catalog# ABX350887-96RXN) and the optical density at 450nm (OD450) was measured by BMG Labtech microplate reader, FLUOstar Omega series.
A third aliquot of quadriceps homogenates was used for global metabolomics analysis using mass spectrometry performed on a Q-TOF Premier mass spectrometer (Waters, Inc.). The metabolite profiling of muscle homogenate was conducted as per the report published previously (3,4). 50 mg of the homogenized tissue was dissolved in chilled extraction buffer (50% methanol in water). The protein was precipitated by addition of equal volumes of acetonitrile; samples were centrifuged. 5 ul of the sample was injected onto a reverse-phase 50 × 2.1mm Acquity 1.7-μm C18 column (Waters Corp, Milford, MA) using an Acquity UPLC system (Waters) with a gradient mobile phase consisting of 2% acetonitrile in water containing 0.1% formic acid (Solvent A) and 2% water in acetonitrile containing 0.1% formic acid (Solvent B) and resolved for 10 min at a flow rate of 0.5 mL/min. The gradient consisted of 100% A for 0.5 min with a ramp of curve 6-100% B from 0.5 to 10 min. The column eluent was introduced directly into the mass spectrometer by electrospray. Mass spectrometry was performed on a Q-TOF Premier mass spectrometer (Waters) operating in either positive-ion or negative-ion electrospray ionization (ESI+/−) mode with a capillary voltage of 3200 V and a sampling cone voltage of 20 V in negative mode and 35 V in positive mode. The desolvation gas flow was set to 800 L/h and the temperature was set to 350 °C. The cone gas flow was 25 L/h, and the source temperature was 120 °C. Accurate mass was maintained by introduction of LockSpray interface of sulfadimethoxine (311.0814 [M+H] + or 309.0658 [M−H]−) at a concentration of 250 pg/μL in 50% aqueous acetonitrile and a rate of 150 μL/min. Data were acquired in centroid mode from 50 to 850m/z in MS scanning. Centroided and integrated mass spectrometry data obtained in duplicates from the UPLC-TOFMS were processed using CEU Mass Mediator database (www.ceumass.eps.uspceu.es).

Result
Quadriceps tissue specific sequencing data was able to map onto the whole mouse genome at the average rate of 20,000 reads per gene. As a result, we captured 10,874 genes (~54% of whole mouse genes) with reads more than 10 counts (counts per million > 0. 5 Functional analysis of DEGs and DECs were seeded to identify 23 significantly enriched (scored by hypergeometric p values) and significantly regulated (measured by Z-score) canonical networks. Table S1 listed all significantly regulated canonical networks. Taking lead from the canonical network analysis, we constructed 6 non-canonical networks, which are enlisted in Table S2.
Upstream regulators defined as any molecular entity that is known to affect the expression of a single gene or a family of genes via direct or indirect causal association (5). Upstream Regulator Analysis (URA) in the Ingenuity Pathway Analysis (IPA) platform was adapted from the Connectivity Map tool (6), where the relationships among the genes were mapped using database developed by literature curation. These upstream regulators could have significant clinical potential, since a large number of genes and their biofunctions can be systematically controlled via targeting single upstream regulator. Table S3 listed the genes regulated by the upstream regulators of interest.
Figure legends: Figure S1. Quadriceps samples' quality control check. Post TMM normalization, the normalized counts of individual samples were plotted in the box and whisker plot. Here the box covers the interquartile range (from Q1 to Q3), the middle line across the box represents the mean value and the two ends of the whisker represents the range (from maximum to minimum). The plot shows that TMM normalization achieved a significant homology across the samples from FLT and CTR.
Table legends: Table S1. The list of canonical pathways significantly perturbed by differentially expressed genes and metabolites. The list is sorted based on the individual pathway's z-score. Pathways that scored a positive z-score above 1.5 were considered activated; while those pathways, which scored a negative z-score lower than -1.5 were considered inhibited. There are two activated and twenty one inhibited pathways listed in this table. The log transformed p-values, number of molecules enriching particular pathway and the molecular IDs (gene symbol and metabolite name) are listed as well.