Identi ﬁ cation of differentially expressed genes in mouse paraspinal muscle in response to microgravity

Lower back pain (LBP) is the primary reason leading to dyskinesia in patients, which can be experienced by people of all ages. Increasing evidence have revealed that paraspinal muscle (PSM) degeneration (PSMD) is a causative contributor to LBP. Current research revealed that fatty in ﬁ ltration, tissue ﬁ brosis, and muscle atrophy are the characteristic pathological alterations of PSMD, and muscle atrophy is associated with abnormally elevated oxidative stress, reactive oxygen species (ROS) and in ﬂ ammation. Interestingly, microgravity can induce PSMD and LBP. However, studies on the molecular mechanism of microgravity in the induction of PSMD are strongly limited. This study identi ﬁ ed 23 differentially expressed genes (DEGs) in the PSM (longissimus dorsi) of mice which were ﬂ own aboard the Bion M1 biosatellite in microgravity by bioinformatics analysis. Then, we performed protein – protein interaction, Gene Ontology function, and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis for the DEGs. We found that Il6ra, Tnfaip2, Myo5a, Sesn1, Lcn2, Lrg1, and Pik3r1 were in ﬂ ammatory genes; Fbox32, Cdkn1a, Sesn1, and Mafb were associated with muscle atrophy; Cdkn1a, Sesn1, Lcn2, and Net1 were associated with ROS; and Sesn1 and Net1 were linked to oxidative stress. Furthermore, Lcn2, Fbxo32, Cdkn1a, Pik3r1, Sesn1, Net1, Il6ra, Myo5a, Lrg1, and Pfkfb3 were remarkably upregulated, whereas Tnfaip2 and Mafb were remarkably downregulated in PSMD, suggesting that they might play a signi ﬁ cant role in regulating the occurrence and development of PSMD. These ﬁ ndings provide theoretical basis and therapeutic targets for the treatment of PSMD.


Introduction
Lower back pain (LBP), including nociceptive, neuropathic, nociplastic, or non-specific pain, is one of the most common musculoskeletal disorders and might be experienced by people of all ages (1-3). As a global public health problem, LBP is known to impose a huge socioeconomic burden worldwide each year (3,4). The maintenance of spinal stability depends on the coordinated movement of the active subsystem [paraspinal muscles (PSM) and tendons], the passive subsystem (vertebrae, intervertebral discs, and ligaments), and the nervous system (5,6). The structure and the function disorder of any of the above components might contribute to LBP. The psoas major, multifidus, and erector spinae are the most important muscle groups of the PSM, which exert a critical role in maintaining the upright posture and the sagittal balance of the spine (5,7,8). In recent years, increasing evidence have revealed that paraspinal muscle (PSM) degeneration (PSMD) is closely correlated with LBP and spine degenerative diseases, such as intervertebral disc degeneration and scoliosis (9-13). Microgravity is caused by factors such as the residual atmosphere in space, which is reported to affect the metabolism of human muscles and bones by regulating the gene expression and the molecular signaling pathways. The space environment is a microgravity environment; this environment can induce PSMD and intervertebral disc degeneration. Thus, microgravity is also a key contributor to LBP (14)(15)(16). Nevertheless, the role and the molecular regulatory mechanisms of microgravity in the induction of PSMD remain largely unknown.
Up to now, most studies focused on multifidus degeneration and spine pathophysiology (9, 10,[17][18][19][20]. Notably, the degeneration of the erector spinae usually occurs earlier, and the degree of degeneration is more serious than that of the multifidus (5,21). The erector spinae consists of the iliocostalis lumborum, longissimus dorsi, and spinous muscles. The developing fatty infiltration, tissue fibrosis, and muscle atrophy have long been considered as characteristic pathological alterations of PSMD (9). Muscle atrophy is associated with abnormally elevated oxidative stress, reactive oxygen species (ROS), and inflammation (15,22). Altered gene expression appears to be a hallmark pathological feature of PSMD. Yang et al. (23) found that the serum CXC chemokine ligand 10 concentrations were significantly increased in patients with PSMD compared with patients without PSMD. Kudo et al. (19) observed that anti-muscular dystrophy gene peroxisome proliferator-activated receptor gamma coactivator 1 alpha (PGC-1a) as well as pro-inflammatory cytokines TNF-a and IL-6 were significantly increased in multifidus from patients with lumbar kyphosis (a reduced lumbar lordosis) compared with normal lumbar lordosis. The genes of the PSM were also differentially expressed in microgravity. Yamakuchi et al. (24) found that 42 genes were differentially expressed in the PSM of rats after 14 days of space flight, among which the heat shock protein 70 and t complex polypeptide 1 were increased, whereas myocyte-specific enhancer binding factor 2C, aldolase A, and muscle ankyrin were decreased. Mirzoev et al. (25) found that ubiquitin ligase MURF-1 was upregulated in the longissimus dorsi of mice after a flight of 30 days. Ogneva et al. (26) found that the content of a-actinin-1/4 and b-actin, respectively, were decreased in the skeletal muscle of mice after a space flight on board the BION-M1 biosatellite for 30 days. However, the effects of microgravity on the gene expression of the longissimus dorsi (erector spinae) remain unclear.
Gambler and colleagues studied the global gene expression profile in the longissimus dorsi of C57BL/N6 male mice after a space flight of 30 days and reported that several genes were associated with insulin resistance and PSM metabolism (27). In this study, we re-analyzed this microarray dataset (GSE94381) downloaded from the public open Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo) (28). R software was used to identify the differentially expressed genes (DEGs). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis were performed to predict the potential biological functions of the DEGs.

Analysis of the mRNA microarray dataset
The GSE94381 microarray dataset was downloaded from GEO database through R software GEOquery package (29). The probes corresponding to multiple molecules for one probe were removed, and we kept only the probes with the largest signal value when encountering probes corresponding to the same molecule. The data was then standardized using the normalizeBetweenArrays function of the limma package in R software. Limma package was also used to identify the DEGs according to the following criteria: -log 10 (P-value) > log 10 20 and |log 2 fold-change (FC)| > 1. The Ggplot2 package in R software was used to analyze and visualize the clustering situation between the different groups (30). The Ggplot2 package was also used to visualize the principal component analysis (PCA), uniform manifold approximation and projection (UMAP), and volcano plots. The heat map was visualized using the complexheatmap package (31).

Venn diagram analysis
In the GSE94381 microarray dataset, the mice were divided into three groups: Bion-flown (BF), Bion ground (BG), and flight control (FC) group (28). The DEGs between the BG and BF groups, the FC and BF groups, as well as the FC and BG groups were intersected to select overlapping genes using http://www. bioinformatics.com.cn, a free online platform for data analysis.

GO and KEGG pathway enrichment analysis
The potential biological functions of the DEGs were predicted through GO and KEGG enrichment analysis. The GO database divides the gene functions into three parts: cellular component (CC), molecular function (MF), and biological process (BP). CC is used to describe the location of gene products in cells, such as endoplasmic reticulum or nucleus; MF mostly refers to the functions of individual gene products, such as binding activity or catalytic activity; while BP mostly refers to an orderly biological process with multiple steps, such as cell growth, apoptosis, and signal transduction. The KEGG pathway analysis is a comprehensive database that integrates genomic, chemical, and systemic functional information. KEGG has four major categories and 17 sub-databases, one of which is called the KEGG pathway. The results of the GO and KEGG analyses were visualized using an online tool (http://www. bioinformatics.com.cn).

Protein-protein interaction network construction
The Search Tool for the Retrieval of Interacting Genes (STRING) database (https://cn.string-db.org/) was widely utilized to assess protein-protein interactions (PPIs) in functional protein association networks (34). Then, we input the DEGs into the multiple protein section of the STRING database and set the organisms as "Mus musculus" to construct the PPI network. The "required score" was set at >0.4. Finally, we generated and downloaded the results of the PPI network analysis.

Statistical analysis
The expression data of several key genes in the GSE94381 dataset were analyzed and visualized using GraphPad Prism software. The statistical significance between the two groups was compared by an unpaired Student's t-test, whereas the differences among more than two groups were evaluated by one-way analysis of variance followed by Turkey's multiplecomparisons test. The results were expressed as mean ± standard deviation. P-value <0.05 was determined to statistical significance (** represents p < 0.01 and *** represents p < 0.001).

Evaluation of the reasonableness of the sample and grouping
To explore the DEGs that can regulate PSMD in microgravity, we analyzed the microarray dataset GSE94381. The basic information of GSE94381 is shown in Table 1. Given that the BF group of mice (experimental group) was flown aboard the Bion M1 biosatellite in microgravity, the authors designed the BG group wherein mice were housed in the same condition but exposed to Earth's gravity and the FC group wherein mice were housed in a standard animal facility to rule out the effect of housing conditions. We analyzed the transcriptome of the BG-BF group, the FC-BF group, and the FC-BG group, respectively.
PCA is a technique for simplifying datasets, which reflects the difference and the distance between different samples by presenting multiple sets of data on the coordinate axis. The closer the distance in the PCA diagram, the more similar the sample composition. As shown in Figures 1A-C, the distance between the same group of samples was closer, while the samples between the BG and BF groups, the FC and BF groups, as well as the FC and BG groups were separated, suggesting that the within-group differences were smaller, while the differences between groups were obvious, and there may be more meaningful results in the subsequent difference analysis. The UMAP plot could clearly distinguish the different groups, which further supported the rationality of sample grouping ( Figures 1D-F). The box plot was utilized to confirm the distribution trends of the hybridization data and the degree of dispersion. As shown in Figures 1G-I, we did not observe

Identification of the DEGs
In our study, the thresholds of -log 10 (P-value) > log 10 20 and | log2 (FC)| > 1 were used to select the DEGs. The volcano plots were used to evaluate the gene expression variation among different groups. The cluster heat map showed the first 20 DEGs. A total of 27 DEGs were downregulated and 40 DEGs were upregulated between the BG and BF groups (Figures 2A, D), 16 DEGs were downregulated and 30 DEGs were upregulated between the FC and BF groups ( Figures 2B, E), and 11 DEGs were downregulated and 21 DEGs were upregulated between the FC and BG groups ( Figures 2C, F).
First, a total of 32 DEGs were identified between two control groups. Second, we observed that there were only five overlapping DEGs by intersecting BG-BF and FC-BG as well as 10 DEGs by intersecting FC-BF and FC-BG through Venn analysis ( Figure 3A), suggesting that the BION-M1 biosatellite, on its own, hardly affected the alteration of gene expression in flown mice muscle compared with the ground control (BG). We also found 23 overlapping DEGs by intersecting BG-BF and FC-BF ( Figure 3A); then, we determined to study the functions of these DEGs. pathway enrichment analysis showed that the 23 DEGs were enriched in the HIF-1 signaling pathway, FoxO signaling pathway, bladder cancer, JAK-STAT signaling pathway, endometrial cancer, renal cell carcinoma, p53 signaling pathway, melanoma, non-small cell lung cancer, and glioma ( Figure 3B). Furthermore, Il6ra, pik3r1, and Cdkn1a were involved in the regulation of JAK-STAT signaling pathway and HIF-1 signaling pathway, thereby mediating cell cycle, proliferation, apoptosis, inflammatory response, etc. (Figures 3C, D, F). Sesn1 and Cdkn1a might regulate DNA repair and damage prevention and cell cycle arrest through mediating the p53 signaling pathway (Figures 3C, E). Additionally, Fbox32, pik3r1, and Cdkn1a were involved in the regulation of FoxO signaling pathway, which can regulate muscle atrophy, cell cycle, etc. (Figures 3C, G).

GO functional enrichment analysis
To further predict the functions of DEGs, we conducted GO functional enrichment analysis. GO function annotations included BP, MF, and CC. As shown in Figures 4A, B, the top 10 BP terms were enriched in the negative regulation of osteoclast differentiation, cellular response to radiation, regulation of smooth muscle cell proliferation, response to reactive oxygen species, negative regulation of myeloid leukocyte differentiation, response to steroid hormone, etc., of which lipocalin-2 (Lcn2), sestrin 1 (Sesn1), and neuroepithelial cell transforming 1 (Net1) were involved in the regulation of cell response to reactive oxygen species, interleukin 6 receptor alpha (Il6ra), phosphatidylinositol 3-kinase, regulatory subunit polypeptide 1 (Pik3r1), and cyclindependent kinase inhibitor p21 (Cdkn1a) were involved in the Frontiers in Endocrinology frontiersin.org groups, of which sesn1 was related to "GATOR2 complex, TORC2 complex, TOR complex, and Seh1-associated complex", and Pik3r1 and Cdkn1a were related to "transferase complex, transferring phosphorus-containing groups" (Figures 4C, D). The genes involved in MF are as follows: calmodulin binding, cytokine receptor binding, SNARE binding, growth factor receptor binding, 1-phosphatidylinositol-3-kinase activity, carbohydrate phosphatase activity, insulin receptor substrate binding, and cyclin-dependent protein serine/threonine kinase inhibitor activity-of which Il6ra was involved in cytokine receptor binding and growth factor receptor binding, and Pik3r1 was related to "calmodulin binding, cytokine receptor binding, growth factor receptor binding, 1-phosphatidylinositol-3-kinase activity, and insulin receptor substrate binding" (Figures 4E, F).

Identification of PSMD-related DEGs
Given that the abnormally increased oxidative stress, ROS, and inflammation and muscle atrophy are the important pathological mechanism of PSMD, the aim of this study was to explore the relevant genes. Ferroptosis is iron-dependent cell death, and the production of a large number of ROS is its hallmark pathological feature (35). A total of 567 FRGs were identified from the FerrDb database via intersecting ferroptosis Frontiers in Endocrinology frontiersin.org driver, suppressor, marker, and unclassified regulator ( Figure 5A). Thus, we merged FRGs, OSRGs, and DEGs to determine the PSMD-related DEGs. The result unveiled that Cdkn1a and Lcn2 were related to ROS and that Sesn1 and Net1 were associated with oxidative stress as well ( Figure 5B). Using the STRING database, we constructed a PPI network of 23 DEGs and visualized them ( Figure 5C). Combining literature reports and the results of functional enrichment analysis, we determined that Il6ra, TNF-a-inducible protein 2 (Tnfaip2), myosin Va (Myo5a), Sesn1, Lcn2, leucine-rich a-2 glycoprotein1 (Lrg1), and Pik3r1 were inflammatory genes; Fbox32, Cdkn1a, Sesn1, and musculoaponeurotic fibrosarcoma oncogene (MAF)/MAF family B (Mafb) were associated with muscle atrophy; Cdkn1a, Sesn1, Lcn2, and Net1 were associated with ROS; and Sesn1 and Net1 were linked to oxidative stress ( Figure 5C).

Discussion
Spinal degeneration and instability are not only due to changes in bone and intervertebral disc structure and function but also to PSMD. The limitation of the current research on the PSM is that they only measure the morphological parameters of PSM on T2-weighted magnetic resonance imaging, including cross-sectional area, fat infiltration rate, etc. The understanding of the molecular pathological mechanism of PSMD is very limited. This year, Anderson et al. (36) collected multifidus muscle for biopsy and confirmed that pro-fibrogenic and proatrophic genes were upregulated and that anti-fibrogenic and This study re-analyzed the gene expression profile of the longissimus dorsi of mice exposed to microgravity for 30 days and found that several DEGs were correlated to oxidative stress, ROS, inflammation, and muscle atrophy, which are the pathogenesis of PSMD. PSM atrophy is a health problem that is currently receiving increasing attention. Identifying the relevant molecular mechanisms and blocking or alleviating PSM atrophy is therefore of utmost importance, yet the mechanisms leading to muscle atrophy are largely unclear. PSM atrophy may be driven by changes in molecules linked to oxidative stress, ROS, and inflammation (9, 15,22,37). Growing evidence has revealed that microgravity can alter the expression levels of muscle atrophyrelated genes (15,(24)(25)(26). To further study the molecular mechanisms of muscle atrophy in response to microgravity, we performed bioinformatics analysis. Muscle-specific ubiquitin E3 ligase Fbox32 (atrogin-1) was reported to be highly expressed in muscle atrophy and inflammation (38). Anderson et al. (36) demonstrated that FBXO32 was found to be remarkably increased in multifidus muscle from 59 patients with lumbar spine pathology by quantitative polymerase chain reaction. Our study also predicted that Fbox32 was upregulated under microgravity, which was consistent with Gambler's result (27). Another E3 ligase muscle atrophy F-Box (MAFb) was downregulated under microgravity (22,27). Mahmassani et al. (39) demonstrated that the expression of MAFB was downregulated in older adults and associated with the change in leg muscle mass. These findings are consistent with the results of our analysis. Additionally, the cyclin-dependent kinase inhibitor P21 (Cdkn1a) is an inducer of cell cycle arrest and senescence-associated gene, which also serves as a mediator in controlling muscle atrophy (40,41) and inflammation (42). Kumar et al. (37) demonstrated that microgravity can substantially induce Cdkn1a expression, which is in agreement with our data. Sestrin 1 (Sesn1) was verified to be significantly downregulated in several models of muscle atrophy, including sarcopenia, and protect muscles against aging-induced atrophy (43). However, Sesn1 was significantly upregulated in mouse PSM in response to microgravity, but its roles and mechanism deserve further study. Thus, Sesn1, Cdkn1a, Fbox32, and Mafb might play a pivotal role in muscle atrophy under microgravity. Chronic muscle inflammation may contribute to the rapid loss of muscle mass, function, and myofibrillar proteins. Our study observed that Il6ra, Pik3r1, and Lrg1 were involved in cytokine receptor binding; Cdkn1a, Il6ra, Pik3r1, and Pfkfb3 were found as well to belong to the HIF-1 signaling pathway via GO and KEGG pathway enrichment analysis. Furthermore, Li et al. (44) reported that the inflammatory gene PIK3R1 was remarkably overexpressed in the total T cells in COVID-19 patients. Liu et al. (45) demonstrated that leucine-rich a-2 glycoprotein 1 (LRG1) may prevent renal fibrosis by inhibiting the inflammatory response and pro-fibrotic cytokines. Interleukin (IL)-6 is a common pro-inflammatory cytokine that regulates many inflammatory pathways. IL6RA is an IL6 receptor. Frempah et al. (46) demonstrated that the absence of IL6RA in keratinocytes can promote an inflammatory response. Wang et al. (47) confirmed that the inhibition of PFKFB3 significantly attenuated the inflammatory response in human valve interstitial cell. Faust et al. (42) found that the decrease in Cdkn1a expression paralleled the reduction in cartilage tissue inflammation, suggesting that Cdkn1a may be associated with cartilage inflammation in osteoarthritis. TNFAIP2, a TNF-ainduced protein, which was an inflammatory-responsive molecule for TNF-a and was upregulated in spinal cord ischemia/reperfusion injury patients (48). Additionally, the expression of protease-activated receptors (PAR) is associated with inflammation, and MYO5A is a PAR1-dependent transcript, which has been implicated in the inflammatory process (49). LCN2, an inflammation-related factor with elevated expression in the skeletal muscle of ob/ob mice with sarcopenia, is associated with muscle atrophy-related inflammation and oxidative stress (50). Keping et al. (51) demonstrated that SESN1 inhibited oxidized low-density lipoprotein-induced activation of NK-kB signaling and reduced the expression of pro-inflammatory cytokines in macrophages. Collectively, the abovementioned studies unveiled that Pfkfb3, Il6ra, Tnfaip2, Myo5a, Sesn1, Lcn2, Irg1, and Pik3r1 were inflammatory genes, suggesting that they might exert a key role in the pathogenesis of PSMD.
Other causes of inducing muscle atrophy are oxidative stress and ROS, which were confirmed to promote muscle protein catabolism and inhibit protein synthesis signaling pathways (52)(53)(54). In addition to mediating inflammation and muscle atrophy, Cdkn1a was also involved in the regulation of oxidative stress and ROS under microgravity (37). Sesn1 belongs to the family of stress-inducible proteins that regulate oxidative stress and ROS, thus protecting the cells from various stimuli (55). Additionally, Lcn2 was confirmed to promote mitochondrial ROS production and alleviate mitochondrial oxidative phosphorylation in rat primary cardiomyocytes (56). However, whether they mediate PSMD via regulating ROS and oxidative stress is an interesting problem which needs future experiments to be verified.
Our study recognizes some limitations. First, the data was downloaded from GEO database, and we did not perform RNA sequencing. Second, we did not perform molecular biology and animal experiments to demonstrate the expression and the roles of the 12 key DEGs in longissimus dorsi muscle in mice. Lastly, there are inherent differences between mice and humans, so the generalization of mouse findings to humans is limited.

Conclusion
In conclusion, a total of 23 DEGs were observed to be remarkably dysregulated in the PSM of mice in microgravity by bioinformatics analysis. Moreover, we found that Il6ra, Tnfaip2, Myo5a, Sesn1, Lcn2, Lrg1, and Pik3r1 were linked to inflammatory response; Fbox32, Cdkn1a, Sesn1, and Mafb were associated with muscle atrophy; Cdkn1a, Sesn1, Lcn2, and Net1 were associated with ROS; and Sesn1 and Net1 were linked to oxidative stress. Notably, Sesn1 was involved in the regulation of inflammatory response, muscle atrophy, ROS, and oxidative stress. These results suggest that they might be a therapeutic candidate target for PSMD treatment.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.