key: cord-0026924-k2b3sty7
authors: McFarland, Amelia J.; Ray, Pradipta R.; Bhai, Salman; Levine, Benjamin D.; Price, Theodore J.
title: RNA sequencing on muscle biopsy from a 5‐week bed rest study reveals the effect of exercise and potential interactions with dorsal root ganglion neurons
date: 2022-02-08
journal: Physiol Rep
DOI: 10.14814/phy2.15176
sha: f54bd19ba43075b2b6491edcaa81c97d6d57765d
doc_id: 26924
cord_uid: k2b3sty7

Sedentary lifestyle, chronic disease, or microgravity can cause muscle deconditioning that then has an impact on other physiological systems. An example is the nervous system, which is adversely affected by decreased physical activity resulting in increased incidence of neurological problems such as chronic pain. We sought to better understand how this might occur by conducting RNA sequencing experiments on muscle biopsies from human volunteers in a 5‐week bed‐rest study with an exercise intervention arm. We also used a computational method for examining ligand–receptor interactions between muscle and human dorsal root ganglion (DRG) neurons, the latter of which play a key role in nociception and are generators of signals responsible for chronic pain. We identified 1352 differentially expressed genes (DEGs) in bed rest subjects without an exercise intervention but only 132 DEGs in subjects with the intervention. Among 591 upregulated muscle genes in the no intervention arm, 26 of these were ligands that have receptors that are expressed by human DRG neurons. We detected a specific splice variant of one of these ligands, placental growth factor (PGF), in deconditioned muscle that binds to neuropilin 1, a receptor that is highly expressed in DRG neurons and known to promote neuropathic pain. We conclude that exercise intervention protects muscle from deconditioning transcriptomic changes, and prevents changes in the expression of ligands that might sensitize DRG neurons, or act on other cell types throughout the body. Our work creates a set of actionable hypotheses to better understand how deconditioned muscle may influence the function of sensory neurons that innervate the entire body.

Skeletal muscle is highly adaptive and undergoes cellular physiology and phenotype changes as well as transcriptomic reprogramming, in particular in response to extended reductions in physical activity Patterson et al., 2018) . Accordingly, in situations whereby mechanical stimuli and load are significantly reducedsuch as prolonged bed rest or immobilization, exposure to microgravity environments, or a highly sedentary lifestyle -deconditioning responses are observed (Baldwin, 1996; Booth et al., 2017; Brocca et al., 2012) . Inactivity and reduced physical activity have been associated with an array of maladaptive responses, with sedentary time considered an important, independent contributor to the development of metabolic and cardiovascular disease, as well as all-cause mortality (Dempsey et al., 2016; Gabriel & Zierath, 2019; Wilmot et al., 2012) . On the other hand, exercise has clear health benefits, which include cardiovascular and metabolic systems and extend to neurological outcomes, such as pain relief (Sluka et al., 2018) and protection from the development of chronic pain (Grace et al., 2016; Leung et al., 2016) .

Following reductions in loading, skeletal muscle experiences loss of mass and fiber size, muscle fibers shift toward a fast phenotype, myofibril content decreases, reduced synthesis of muscle protein occurs, and muscle protein is degraded by proteasomal and autophagy mechanisms (Adams et al., 2003; DeMartino & Ordway, 1998; Gallagher et al., 2005; Haddad et al., 2005; Hastings et al., 2012; Krainski et al., 2014b; Phillips & McGlory, 2014) . Previous studies have demonstrated that exercise countermeasures, such as rowing ergometry and resistive strength exercises, which were used here, are capable of preserving cardiac and skeletal muscle mass and limit cardiovascular deconditioning associated with prolonged bed rest (Hastings et al., 2012; Krainski et al., 2014b) . Indeed, the integrated cardiovascular benefits associated with rowing ergometry mean that it has been suggested as a strategy to counter severe deconditioning scenarios, such as prolonged spaceflight (Hastings et al., 2012; Krainski et al., 2014b) , or chronic diseases such as the Postural Orthostatic Tachycardia Syndrome (Fu et al., 2010) .

While the vast majority of studies in this field to date have focused on functional and/or structural changes associated with bed rest or limb immobilization, gene expression profiling tools have also permitted exploration of the transcriptional changes which accompany muscle deconditioning, and how exercise may lessen or counteract these (Pillon et al., 2020) . Microarray technology has been used in both acute limb immobilization as well as acute and chronic bed rest studies to explore associated transcriptional changes (Abadi et al., 2009; Chopard et al., 2009; Fernandez-Gonzalo et al., 2020; Lammers et al., 2012; Urso et al., 2006) , while RNAsequencing/RNA-seq, which allows for a deeper investigation of the transcriptome compared to microarray analysis (Stark et al., 2019) , has only been employed in short term bed rest studies (≤14 days; Mahmassani et al., 2019b Mahmassani et al., , 2021 . Despite the link between physical inactivity and onset of musculoskeletal complications and systemic symptoms, knowledge of the mechanisms that underlie this link is not well understood. Moreover, it is not clear whether deconditioned muscle may contribute to disease mechanisms beyond the cardiovascular or musculoskeletal systems, for instance through an action on neurological systems (Booth et al., 2017) , which could contribute to the somatic hypervigilance that is present in deconditioning syndromes (Cutsforth-Gregory & Sandroni, 2019) . We have recently developed computational methods that allow for detailed examination of pharmacological interactions between any tissue and human dorsal root ganglion (DRG) neurons . These neurons play a critical role in pain signaling and tactile and thermal sensation, and signals derived from muscle may interact with these neurons to mediate neurological effects of exercise or deconditioning. Other cell types may be influenced by signals from muscle and the framework we establish in this study could be used to interrogate interactions with any number of cell types ranging from motor neurons to immune cells.

The purpose of the present study was to use RNA sequencing to gain mechanistic insight into transcriptional changes which occur following muscle deconditioning in long-lasting bed rest. We analyzed muscle biopsy samples from healthy individuals which had been previously collected as part of a 5-week head-down-tilt bed rest investigation (Krainski et al., 2014b) . Using a discovery-based approach, we performed bulk RNA-seq on muscle samples and sought to identify genes, pathways and processes which were differentially expressed between individuals who underwent complete bed rest, and those who performed a high-intensity exercise routine during their bed rest. Additionally, we evaluated the potential impact of transcriptional changes in muscle on DRG neurons using a genome-wide ligand-receptor pair database curated for pharmacological interactions relevant to neuro-immune systems . This work represents the longest duration of human physical inactivity analyzed via RNA sequencing technology and integrates this new knowledge into the context of potential impacts on DRG neurons.

Skeletal muscle biopsies were collected at the Institute for Exercise and Environmental Medicine under approved IRB protocol STU 102010-070. Methods and subject characteristics for this study have been previously described (Krainski et al., 2014b) shown in Table 1 and in Supplementary Data, Sheet 1. The samples were collected as part of the (Krainski et al., 2014a) study. Subjects were healthy, nonsmoking, adults aged 20-54 years. Subjects were randomly assigned to one of two groups: bed rest only (BR-CON; n = 9, 1 female), or bed rest with exercise countermeasure (BR-EX; n = 16, 2 female). All subjects underwent −6-degree head down tilt bed rest for 35 consecutive days. BR-EX subjects received an exercise countermeasure consisting of rowing ergometer training on 6 days/week and biweekly resistance training for the duration of the bed rest period. Details of the exercise countermeasure, effects on muscle mass and performance and other metrics are all available (Krainski et al., 2014a) . A needle muscle biopsy of 300 mg tissue on average was obtained from the mid-belly of the right (pre) and left (post) vastus lateralis muscle before and after the 5-week bed rest period (Krainski et al., 2014b) . Part of each biopsy was used for a previously published study (Krainski et al., 2014b) ; at least 50 mg per sample remained for use in the present study. The samples were frozen upon collection and stored at −80°C until used for RNA sequencing for this study. 

Based on available muscle following initial ultrastructural and histochemical analysis (Krainski et al., 2014a) biopsy samples from 13 subjects underwent RNA sequencing, consisting of all three female samples (1 BR-CON, 2 BR-EX), and ten male samples (5 BR-CON, 5 BR-EX; matched for age and exercise tolerance at baseline). Researchers performing sequencing and subsequent analyses were informed of subject pair matches and whether samples were from baseline or post-intervention, but were blinded to intervention received by individual subjects (i.e. the researchers did not know whether subjects were sedentary or exercised during bed rest until the first round of analysis was complete). Following RNA extraction (RNeasy Plus Universal Mini Kit, Qiagen), RNA yield was quantified using a Nanodrop system (ThermoFisher Scientific), and RNA quality was assessed by fragment analyzer (Advanced Analytical Technologies). Stranded mRNA library preparation and sequencing were completed at the University of Texas at Dallas (UT Dallas), TX, USA. 76bp, single-end sequencing of RNA-seq libraries was performed on the Illumina Hi-Seq sequencing platform. The samples were run in a blinded fashion.

RNA-seq read files (fastq files) were checked for quality by FastQC (Babraham Bioinformatics) and read trimming was done based on Phred score and per-base nucleotide content, with trimmed, 60 bp reads from the 13-72 bp positions being used for downstream analysis. Trimmed reads were mapped to the reference transcriptome GENCODE v27 (Frankish et al., 2019) , and its underlying reference genome GRCh38.p10 retaining only uniquely mapped reads by using the STAR (v2.7.3a) aligner (Dobin et al., 2013) with the following command line mapping parameters:

AllBestScore --outBAMcompression -1 --outBAMsortingThreadN 20 --outFilterType BySJout --outFilterMulti-mapNmax 10 --outFilterMismatchNoverReadLmax 0.06 --twopassMode Basic --quant-Mode GeneCounts --alignSJoverhangMin 5 --alignSJDBoverhangMin 3 --alignIntron-Max 1000000 --alignMatesGapMax 1000000

Relative abundance per sample (in Transcripts per Million or TPM) was calculated from the mapped reads using Stringtie (Pertea et al., 2015) . Genes were filtered to only retain protein-coding genes and TPMs were renormalized per sample to a million based only on the set of coding genes A fold-change ratio (post-intervention/ baseline) was calculated for each coding gene in every subject (Supplementary Data, Sheet 2). Data is available for download through GEO, accession number GSE18 6045.

The set of subjects was initially partitioned into pairs of subjects with matching sex and comparable ages (Table  1) . For each subject pair, exactly one was BR-CON, while the other was BR-EX. Blinded researchers were told which sample belonged to which pair, but which member of the pair received which intervention was not revealed to researchers. The female samples were not included as part of the blinded analysis.

Hierarchical clustering of baseline and postintervention samples was performed separately using TPMs of all protein-coding genes, using (1 − Pearson's correlation coefficient) as the distance metric ( Figure  1a ,b). Multiple strategies were utilized to identify potential genes of interest from the blinded analysis of the dataset. First, genes with the potential for differential regulation in the control and exercise groups were identified. This was achieved through identifying genes within each subject pair where one subject had a post-intervention/baseline absolute log 2 fold-change of ≥0.585 (corresponding to a fold change of ≥1.5 or ≤0.66), while the other subject had either no change (fold change between 0.66 and 1.5) or a fold change trend in the opposing direction. Genes that consistently appeared in the list for each subject pair were identified by intersecting the gene lists for all 5 subject pairs (Supplementary Data, Sheet 3). Additionally, genes with the potential for differential expression between baseline and post-intervention samples, irrespective of the intervention, were also identified. Lastly, genes with the potential for dysregulation due to bed rest were identified. This was achieved through calculating the relative variance (ratio of variance to mean) of TPMs values across sample pairs; the top 500 genes with the highest relative variance identified for each sample pair. Genes common to the top 500 highest relative variance genes in all sample pairs were also identified (Supplementary Data, Sheet 4). The relative variance and fold change gene lists generated through the blinded analysis were then compared, and genes that appeared in both lists were used to create a final gene set (Supplementary Data, Sheet 5).

To predict the intervention groupings for the male subjects, hierarchical cluster analysis of the TPMs for the fold-change gene list (Supplementary Data, Sheet 3; Figure 1c ) in baseline and post-intervention samples together were used and were also contrasted with previously described hierarchical cluster analysis of protein-coding gene TPMs post-intervention ( Figure  1b ). Post-intervention samples for four subject pairs (pairs 1-4) were segregated into two clusters (Figure 1c ), suggesting which subjects may have received the same intervention, Post-intervention samples for the fifth pair (pair 5) were outliers in the clustering, and 5A and 5B were manually assigned to the two groups based on gene expression trends. Based on the predicted subject groupings, which group may have received which intervention was also predicted through comparison of this intersected gene set to previous bed rest datasets (Brocca et al., 2012; Mahmassani et al., 2019a) .

Coding of subjects and samples (Table 1 ) was kept at the Levine lab while sequencing analysis was done in the Price lab. The Levine lab revealed the coding only after blinded data analysis from the Price lab was completed.

After completion of the blinded analysis, researchers were unblinded to the interventions of each subject. Differentially expressed genes reported in the formal analysis in this report were identified using the R/Bioconductor package DEseq2, a well-established differential expression analysis tool that models read counts using a negative binomial distribution, and models dispersion as a function of the expression level (Love et al., 2014) . Data were analyzed using a nested design to identify group-specific condition effects for BR-EX and BR-CON subjects post-intervention versus baseline. Pairing was done for within-subject preand post-biopsies to control for volunteer-specific genomic landscape differences as well as medical or other lifestyle differences that cannot be controlled in a small sample size, such as in this study. Protein-coding genes with a TPM > 0.1, and an adjusted p-value of <0.05 were used to identify differentially expressed gene sets for each condition (Supplementary Data, Sheets 6 and 7). p-value adjustment in DEseq2 is based on controlling the false discovery rate to <=0.05 using the Benjamini-Hochberg correction, using the p.adjust function in limma/R (Wettenhall & Smyth, 2004) . Gene set enrichment analysis on the differentially expressed gene (DEG) sets was performed using the Enrichr framework (Kuleshov et al., 2016) . DEGs with potential for effects on DRG neurons were identified using a genome-wide ligand-receptor pair database curated for pharmacological interactions. An interactome F I G U R E 1 Dendrogram based on hierarchical clustering of TPMs from male subjects (n = 10) with predicted groupings in blue and green. (a) All protein-coding genes from paired subjects at baseline. (b) All protein-coding genes from paired subjects post-intervention. (c) Protein-coding genes meeting fold-change criteria (as identified from the blinded analysis) from paired subjects both at baseline and postintervention. The distance (calculated as 1 -Pearson's correlation coefficient between the vectors of coding gene TPMs in two samples) fitted to an ultrametric tree is shown on the x-axis was generated through intersecting the list of DEGs and their expression levels with a ligand-receptor pair list as previously described .

Alternative splicing analysis was only performed for the exon retention/skipping event in PGF that caused gain/ loss of heparin-binding functionality in the encoded peptide respectively. Uniquely mapped bridge reads on the correct strand were counted using IGV (Thorvaldsdottir et al., 2013) for the relevant exon-exon junctions. Reads that overlapped the relevant exon that was identified this way were pooled for control and exercised subjects separately, and normalized to calculate the relative abundance of the alternatively spliced exon in question (in Reads per Kilobase per Million mapped reads). The number of mapped reads that were used for the normalization was the number of unambiguously mapped, genic reads that were mapped by the STAR aligner (Dobin et al., 2013) (relevant transcript and exon ids are provided in Supplementary Data, Sheet 1).

Of 58,763 total genes in the reference genome, we detected 36,072 genes in our samples. All samples yielded high-quality raw sequence data. After filtering for proteincoding genes only, 17,730 unique genes remained.

Hierarchical cluster analysis of all male baseline samples showed high homology across all samples, and between some pairs (Figure 1a) . Post-intervention, several non-paired subjects were seen to cluster together (Figure 1b) , including subjects 4A and 5B, subjects 3B and 4B, subjects 1A and 2B, and subjects 2A and 1B. Through the use of our differential fold-change algorithm, where subjects within a pair demonstrated inconsistent or divergent gene expression post-intervention, we identified 263 genes that were consistently dissimilar between all subject pairs (Supplementary Data, Sheet 3). Additionally, we identified genes that were consistently trending in the same direction after bed rest irrespective of intervention (Supplementary Data, Sheet 3). Calculation of relative TPM variance across all subject pairs identified 141 genes which were consistently in the top 500 highest relative variance genes (Supplementary Data, Sheet 4). There were 17 genes that appeared consistently on all blinded analysis gene lists (Supplementary Data, Sheet 5).

Hierarchical cluster analysis of gene lists generated through fold change and relative variance computations similarly suggested high homology between most samples, particularly paired subjects, at baseline; however, post-intervention samples were observed to cluster primarily in non-paired groupings, consistent with previous hierarchical cluster analysis of all gene data. Based on the consistency of these post-intervention non-pair groupings, it was predicted that subjects 1B, 2A, 3B, 4B, and 5A belonged to one group, while subjects 1A, 2B, 3A, 4A, and 5B belonged to the other group. The pattern of expression for the 12 genes consistently expressed in all gene lists of interest were then reviewed for the two groups. Based on expression trends in previous bed rest studies (in particular, XIRP1, TNFRSF12A, and HSPB1; Brocca et al., 2012; Mahmassani et al., 2019a) , it was hypothesized that the control (bed rest only) group consisted of subjects 1A, 2B, 3A, 4A, and 5B, and that the intervention (bed rest plus exercise) group consisted of subjects 1B, 2A, 3B, 4B, and 5A. Following completion of the blinded analysis, investigators were unblinded to the subject grouping and intervention. The blinded analysis correctly predicted the grouping for all ten subjects, but incorrectly predicted the intervention received.

Following unblinding of subject intervention, DESeq2 was used to determine which genes were differentially expressed between control and intervention groups for the male and female subjects. DESeq2 analysis found 1352 differentially expressed genes (DEGs) in subjects in the BR-CON group (591 up-, 761 downregulated post-intervention vs. baseline, p < 0.05; see Supplementary Data, Sheet 6). In BR-EX subjects only 132 DEGs were identified (69 up-, 63 downregulated post-intervention vs. baseline, p < 0.05; see Supplementary Data, Sheet 7). Of the DEGs identified, 36 genes were consistent between the two groups; all other genes were differentially expressed in the BR-CON group were not significantly altered in subjects who received the exercise countermeasure. Twenty-two DEGs had consistent differential expression patterns in both BR-CON and BR-EX groups (16 upregulated in both groups and six downregulated in both groups), and an additional 14 genes had opposing trends between the groups (i.e., were upregulated in one group and downregulated in the other). Of those genes which were differentially expressed in both BR-CON and BR-EX groups, LBP and LRRC52 were both in the top 10 downregulated genes when ranked by magnitude of fold change (Table 2) , while KIT was found to be a top 10 upregulated gene for both groups.

The potential functional implications of gene expression changes were then explored through performing enrichment analysis on the DEGs from each group. In BR-CON, key upregulated pathways identified by functional enrichment analysis included interferon-mediated signaling pathways, transcriptional regulation pathways, and antigen processing. Key pathways enriched by downregulated DEGs in BR-CON related to mitochondrial respiration and electron transport (Table 3 ). In contrast, pathways enriched by downregulated DEGs in BR-EX included transcriptional regulation pathways, cytokine-mediated signaling, and fatty acid oxidation regulation; pathways enriched by upregulated DEGs included multiple pathways relating to endothelial cell development, proliferation, and regulation (Table 4 ).

To assess the potential impact of transcriptional changes in muscle from bed rest subjects on the peripheral nervous system, we used a genome-wide ligand-receptor pair database curated for pharmacological interactions relevant to neuro-immune systems . This interactome permits analysis of a given gene list, and identifies those genes which are capable of ligand-receptor interactions because the ligands have at least one known receptor expressed by human dorsal root ganglia neurons. Upregulated DEGs from BR-CON and BR-EX were analyzed via the interactome. Of the 591 genes upregulated post-intervention in BR-CON, 26 genes were found to have known receptors in human DRG (Table 5) , including PGF, IL18, CCL2, EFNA1, CCN1, and CXCL2. In the BR-EX group, there were seven upregulated genes post-intervention with known DRG receptors (Table 6 ). We looked more closely at PGF isoforms that encode peptides with heparin-binding

Bed rest plus exercise countermeasure (BR-EX) functionality because these interact with neuropilin-1 and neuropilin-2 (Dewerchin & Carmeliet, 2014) , receptors that are highly expressed by DRG neurons (Moutal et al., 2020; Tavares-Ferreira et al., 2021) . We found that these PGF isoforms were robustly detected suggesting potential to interact with nociceptor expressed neuropilin receptors (Reads per Million mapped reads or RPM > 0.5 in 22 out of 26 samples, Supplementary Data, Sheet 1).

The analysis described above yields several key conclusions. Muscle biopsy RNA sequencing shows that 5 weeks of bed rest causes a large transcriptional change in muscle that is mostly prevented by an exercise intervention. The response in deconditioned muscle mostly involves the downregulation of gene expression for regulators of virtually all aspects the mitochondrial activity, an effect that was completely reversed by exercise intervention. This transcriptional change in sedentary muscle also produces a set of ligands in muscle tissue that may interact with receptors on peripheral sensory neurons that provide sensory and nociceptive information from muscle and/or other tissues that could be exposed to factors released from muscle. This muscle to DRG neuron signaling may play an underappreciated role in the physiological sequalae that emerge from muscle deconditioning (Figure 2) . Our data generate a testable set of hypotheses for mediators that may come from deconditioned muscle to sensitize nociceptors and increase pain susceptibility. The findings presented here provide a new perspective on the transcriptional response to bed rest and how this is altered by an exercise intervention. Consistent with other biochemical and physiological measures, exercise during bed rest minimized changes in the muscle transcriptional response. These changes included normalization of gene expression linked to muscle mitochondrial activity, almost entirely downregulated in the no exercise group, and downregulation of genes that are likely involved in inflammation and signaling to the peripheral nervous system, which were upregulated in the no exercise group. Many of the genes that were upregulated in muscle samples without exercise intervention have been reported in previous studies (Abadi et al., 2009; Chopard et al., 2009; Fernandez-Gonzalo et al., 2020; Lammers et al., 2012; Urso et al., 2006) , but this study used a longer bed rest interval, and RNA sequencing provides a broader picture of how gene expression is altered over this time period. A shortcoming is the use of bulk RNA sequencing rather than a single cell approach that could provide more insight into cellular phenotype (Stark et al., 2019) . We decided on a bulk RNA sequencing approach for several reasons. First, the samples were collected years prior and stored at −80°C so cellular integrity may have been compromised in samples. Second, protocols for nuclei isolation are not optimized for samples that have been frozen for years. Finally, we reasoned that a bulk sequencing approach would give us a fuller picture of transcriptional changes that could be informative for the design of future prospective studies. We focused our ligand-receptor interaction on muscle to sensory neuron communication because of the link between exercise and pain. A consistent observation in clinical and preclinical pain studies is that exercise is one of the most successful interventions for chronic pain (Lesnak & Sluka, 2020; Merkle et al., 2020; Sluka et al., 2018) . One potential explanation for this effect is the production of myokines from exercise that act on sensory neurons to decrease their excitability (Sharma et al., 2010) . Another is that lack of exercise, the primary cause of muscle deconditioning, could produce factors from muscle tissue that promote pain through an action on DRG sensory neurons. Our hypothesis-generating analysis provides a framework for testing this latter idea in mechanistic studies. We found a number of ligands, including chemokines (CXCL2 and CXCL9), interleukins (IL18), and growth factors (PGF) that were strongly upregulated in the BR-CON group that have cognate receptors that are expressed by human sensory neurons, including nociceptors that play a key role in detecting pain and are sensitized in chronic pain states (Figure 2 ). IL-18 has been associated with pain promotion in animal models (Verri et al., 2007) and in cancer pain patients (Heitzer et al., 2012) . IL-18 receptors are expressed in human DRG but IL-18 may also act indirectly through endothelin receptors (Verri et al., 2004) , which are robustly expressed by human nociceptors (Tavares-Ferreira et al., (2) Pre-existing human DRG sequencing data and a ligand-receptor interactome (3) enabled looking at ligand-receptor interactions between deconditioned muscle and DRG neurons. (4) Key muscle ligands with increased expression in deconditioned muscle (in red) matched to receptors that are expressed by human DRG nociceptors (in blue). Made with Biore nder.com 2021). Placental growth factor (PGF) is a vascular endothelial growth factor (VEGF) family growth factor that can act via neuropilin receptor (NRP1 and NRP2 genes) complexes (Dewerchin & Carmeliet, 2014) . NRP1 was recently shown to play an important role in promotion of neuropathic pain in rats (Moutal et al., 2020) and is also highly expressed by human nociceptors (Tavares-Ferreira et al., 2021) , suggesting a potential novel role of PGF in pain signaling via this receptor. Finally, while interferon γ was not directly detected as an upregulated gene, downstream interferon γ signaling was the most enriched GO term detected in the BR-CON group. Muscle deconditioning may promote increased interferon γ expression in other cells types that can then act on human nociceptors which express interferon γ receptors (Tavares-Ferreira et al., 2021) . Importantly, interferon γ has been positively associated with several clinical pain conditions (Kamieniak et al., 2020; Parkitny et al., 2013) . These ligand-receptor interactome findings lay a foundation for future hypothesis-based work using human DRG neurons (Renthal et al., 2021) to better understand how deconditioned muscle may signal to the peripheral nervous system.

The ligands that were increased in deconditioned muscle may act on other cell types, in addition to DRG neurons, which were the focus of our study. The increasing availability of RNA sequencing datasets for human organs (Mele et al., 2015) and single cells (Consortium & Quake, 2021) can allow other investigators to use our dataset and interactome framework  or other frameworks; Efremova et al., 2020) to interrogate ligandreceptor interactions between deconditioned muscle and other cell types. As an example, we looked at two of the ligand-receptor interactions noted in our deconditioned muscle -DRG neuron interactome using data from the GTEx Consortium (Mele et al., 2015) , the Tabula Sapiens Consortium (Consortium & Quake, 2021), and a recent preprint with single nucleus RNA sequencing of human spinal cord, including motor neurons (Zhang et al., 2021) , and a complementary spatial RNA sequencing dataset on human spinal cord (Maniatis et al., 2019) . For PGF-NRP1, NRP1 is widely expressed in peripheral tissues (Mele et al., 2015) , and among immune cells is enriched in macrophages and mesenchymal stem cells (Consortium & Quake, 2021) . The receptor is also expressed by motor neurons (Maniatis et al., 2019; Zhang et al., 2021) , therefore this PGF-NRP1 interaction may influence motor neuron, macrophage, and mesenchymal stem cell function, in addition to DRG neurons. For IL18-IL18R1, IL18R1 is also widely expressed in human organs (Mele et al., 2015) , but enriched in a different immune cell population, neutrophils (Consortium & Quake, 2021) . The receptor mRNA is not detected in human motor neurons (Maniatis et al., 2019; Zhang et al., 2021) . As tools and datasets continue to proliferate to enable these types of analyses, it will be interesting to generate more comprehensive views of how ligands from deconditioned muscle may influence broad physiological systems rather than focusing on individual cell types.

There are several additional limitations to the study that deserve mention. First, we have not been able to validate the release of factors from muscle cells using blood or biopsy samples. Our work here sets the stage for such studies in the future. Second, we have focused on effects on peripheral neurons, but some of these factors could potentially also have an effect on CNS neurons if they are able to cross the blood-brain barrier. Another possibility is interactions with peripheral immune cells which could then have an effect on either the peripheral or central nervous system (Grace et al., 2021) . This can be an area of future focus, where immune cell samples could potentially be isolated from the same subjects having muscle biopsy. Of course, we can also not discount the important effect that exercise has on the brain and how this may also influence pain outcomes (Lesnak & Sluka, 2020) . Third, we did not do an interactome analysis on downregulated genes. We did not do this because most of these genes were involved in intracellular mitochondrial signaling and because the impact of downregulated ligands on receptors is not immediately obvious. We note that we did not observe downregulation of pain resolution factors like IL-10 or IL-4 in our deconditioned muscle samples. Downregulation of these ligands could remove endogenous tone that would be expected to inhibit nociceptor sensitization (Durante et al., 2021; Eijkelkamp et al., 2016; Milligan et al., 2005) . Finally, we did not validate results with an alternative method like quantitative polymerase chain reaction (qPCR); however, the utility of such validation methods for RNA-sequencing experiments such as those done here have been called into question (Coenye, 2021) .

The data and analyses described here are a useful resource for understanding how muscle responds to longlasting inactivity. As noted in our results, we were able to accurately classify groups based on their transcriptional response, but we incorrectly called the exercise intervention group prior to unblinding. This mistake was made based on our a priori assumption that antiinflammatory myokines would be produced by exercised muscle. Instead, we discovered that deconditioned muscle produces a broad set of secreted molecules, including a mix of inflammatory and anti-inflammatory mediators, which are likely to have a broad impact on many cell types in the body. Therefore, our unbiased approach reveals new insight into the potential impact that deconditioned muscle may have on the physiology of many organ systems. While our analysis has focused on DRG neurons, our data can be mined for additional contexts. Our results with a relatively small sample size suggest that future prospective studies can utilize single-cell approaches to build substantially on the findings described here.

Limb immobilization induces a coordinate down-regulation of mitochondrial and other metabolic pathways in men and women

Skeletal muscle unweighting: Spaceflight and ground-based models

Effect of spaceflight on the functional, biochemical, and metabolic properties of skeletal muscle

Lack of exercise is a major cause of chronic diseases

Role of inactivity in chronic diseases: evolutionary insight and pathophysiological mechanisms

The time course of the adaptations of human muscle proteome to bed rest and the underlying mechanisms

Largescale mRNA analysis of female skeletal muscles during 60 days of bed rest with and without exercise or dietary protein supplementation as countermeasures

Do results obtained with RNA-sequencing require independent verification? Biofilm, 3, 100043

The Tabula Sapiens: A multiple organ single cell transcriptomic atlas of humans. bioRxiv

Clinical neurophysiology of postural tachycardia syndrome

Ubiquitin-proteasome pathway of intracellular protein degradation: Implications for muscle atrophy during unloading

Sitting less and moving more: improved glycaemic control for type 2 diabetes prevention and management

Placental growth factor in cancer

STAR: Ultrafast universal RNA-seq aligner

Adenosine A3 agonists reverse neuropathic pain via T cell-mediated production of IL-10

Cell PhoneDB: Inferring cell-cell communication from combined expression of multi-subunit ligand-receptor complexes

IL4-10 fusion protein is a novel drug to treat persistent inflammatory pain

Three months of bed rest induce a residual transcriptomic signature resilient to resistance exercise countermeasures

GENCODE reference annotation for the human and mouse genomes

Cardiac origins of the postural orthostatic tachycardia syndrome

Circadian rhythms and exercise-Re-setting the clock in metabolic disease

Effects of 84-days of bedrest and resistance training on single muscle fibre myosin heavy chain distribution in human vastus lateralis and soleus muscles

Prior voluntary wheel running attenuates neuropathic pain

The neuroimmunology of chronic pain: From rodents to humans

Pretranslational markers of contractile protein expression in human skeletal muscle: Effect of limb unloading plus resistance exercise

Effect of rowing ergometry and oral volume loading on cardiovascular structure and function during bed rest

IL-7, IL-18, MCP-1, MIP1-beta, and OPG as biomarkers for pain treatment response in patients with cancer

IFN-gamma correlations with pain assessment, radiological findings, and clinical intercourse in patient after lumbar microdiscectomy: Preliminary study

The effect of rowing ergometry and resistive exercise on skeletal muscle structure and function during bed rest

The effect of rowing ergometry and resistive exercise on skeletal muscle structure and function during bed rest

Enrichr: A comprehensive gene set enrichment analysis web server 2016 update

Expression of genes involved in fatty acid transport and insulin signaling is altered by physical inactivity and exercise training in human skeletal muscle

Mechanism of exercise-induced analgesia: What we can learn from physically active animals

Regular physical activity prevents chronic pain by altering resident muscle macrophage phenotype and increasing interleukin-10 in mice

Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2

Reduced physical activity alters the leucinestimulated translatome in aged skeletal muscle

Age-dependent skeletal muscle transcriptome response to bed rest-induced atrophy

Disuse-induced insulin resistance susceptibility coincides with a dysregulated skeletal muscle metabolic transcriptome

Spatiotemporal dynamics of molecular pathology in amyotrophic lateral sclerosis

Human genomics. The human transcriptome across tissues and individuals

The interaction between pain and movement

Controlling pathological pain by adenovirally driven spinal production of the anti-inflammatory cytokine, interleukin-10

SARS-CoV-2 Spike protein co-opts VEGF-A/Neuropilin-1 receptor signaling to induce analgesia

Inflammation in complex regional pain syndrome: A systematic review and meta-analysis

Sedentary behaviour and risk of all-cause, cardiovascular and cancer mortality, and incident type 2 diabetes: A systematic review and dose response meta-analysis

StringTie enables improved reconstruction of a transcriptome from RNA-seq reads

CrossTalk proposal: The dominant mechanism causing disuse muscle atrophy is decreased protein synthesis

Transcriptomic profiling of skeletal muscle adaptations to exercise and inactivity

Human cells and networks of pain: Transforming pain target identification and therapeutic development

Aerobic exercise alters analgesia and neurotrophin-3 synthesis in an animal model of chronic widespread pain

Exercise-induced pain and analgesia? Underlying mechanisms and clinical translation

RNA sequencing: The teenage years

Spatial transcriptomics reveals unique molecular fingerprints of human nociceptors. bioRxiv

Integrative Genomics Viewer (IGV): High-performance genomics data visualization and exploration

Analysis of human skeletal muscle after 48 h immobilization reveals alterations in mRNA and protein for extracellular matrix components

Antigen-induced inflammatory mechanical hypernociception in mice is mediated by IL-18

Interleukin-18 induces mechanical hypernociception in rats via endothelin acting on ETB receptors in a morphine-sensitive manner

A ligandreceptor interactome platform for discovery of pain mechanisms and therapeutic targets

limmaGUI: A graphical user interface for linear modeling of microarray data

Sedentary time in adults and the association with diabetes, cardiovascular disease and death: Systematic review and metaanalysis

Single-nucleus transcriptomic atlas of spinal cord neuron in human