Prolactin-Releasing Peptide Differentially Regulates Gene Transcriptomic Profiles in Mouse Bone Marrow-Derived Macrophages

Prolactin-releasing Peptide (PrRP) is a neuropeptide whose receptor is GPR10. Recently, the regulatory role of PrRP in the neuroendocrine field has attracted increasing attention. However, the influence of PrRP on macrophages, the critical housekeeper in the neuroendocrine field, has not yet been fully elucidated. Here, we investigated the effect of PrRP on the transcriptome of mouse bone marrow-derived macrophages (BMDMs) with RNA sequencing, bioinformatics, and molecular simulation. BMDMs were exposed to PrRP (18 h) and were subjected to RNA sequencing. Differentially expressed genes (DEGs) were acquired, followed by GO, KEGG, and PPI analysis. Eight qPCR-validated DEGs were chosen as hub genes. Next, the three-dimensional structures of the proteins encoded by these hub genes were modeled by Rosetta and Modeller, followed by molecular dynamics simulation by the Gromacs program. Finally, the binding modes between PrRP and hub proteins were investigated with the Rosetta program. PrRP showed no noticeable effect on the morphology of macrophages. A total of 410 DEGs were acquired, and PrRP regulated multiple BMDM-mediated functional pathways. Besides, the possible docking modes between PrRP and hub proteins were investigated. Moreover, GPR10 was expressed on the cell membrane of BMDMs, which increased after PrRP exposure. Collectively, PrRP significantly changed the transcriptome profile of BMDMs, implying that PrRP may be involved in various physiological activities mastered by macrophages.

The original biological function of PrRP is to promote the release of prolactin from cultured pituitary cells [2]. Then, follow-up research results show that PrRP has an important role in the neuroendocrine system, including reducing food intake in fasting and freeeating rat models [8], regulating energy homeostasis and increasing energy expenditure [9], mediating the anorexia effect of nerves [10], protecting neuron on several rodent neurodegenerative disorder models [11], and participating in stress response [12]. However, the effect of PrRP on macrophages has not been reported yet.

PrRP Demonstrated No Significant Effect on the Shape of BMDMs
As displayed in Figure 2A, the cell morphology of BMDMs was detected by electron microscopy. In addition, the purity of BMDMs was evaluated using flow cytometry with anti-CD11b and anti-F4/80 ( Figure 2B). The results of flow cytometry showed that the double-positive rate of the BMDMs control group (no antibody was added) was 0.14% (the left side of Figure 2B), whereas the double-positive rate of BMDMs (incubated with anti-F4/80 and anti-CD11b) 96.7% (right side of Figure 2B). Hence, these data hinted that the purity of BMDMs was acceptable.

PrRP Demonstrated No Significant Effect on the Shape of BMDMs
As displayed in Figure 2A, the cell morphology of BMDMs was detected by electron microscopy. In addition, the purity of BMDMs was evaluated using flow cytometry with anti-CD11b and anti-F4/80 ( Figure 2B). The results of flow cytometry showed that the double-positive rate of the BMDMs control group (no antibody was added) was 0.14% (the left side of Figure 2B), whereas the double-positive rate of BMDMs (incubated with anti-F4/80 and anti-CD11b) 96.7% (right side of Figure 2B). Hence, these data hinted that the purity of BMDMs was acceptable.
Next, the morphological characters of BMDMs before and after PrRP exposure were examined with an optical microscope. As shown in Figure 2C, BMDMs were oval before PrRP (1 nM) exposure, and the shape of BMDMs did not change remarkably after PrRP exposure for 18 h. Overall, PrRP demonstrated no significant effect on the shape of BMDMs.

Identification of DEGs
To explore the effect of PrRP on the transcriptome of BMDMs, cell samples were detected by RNA-seq sequencing ( Figure 3A-D, Table S1, and Supplementary File 2). The quality control test results showed that our RNA-seq was qualified (Supplementary File S2 Figure S1A and S1B). A total of 410 DEGs were obtained, of which 371 genes were upregulated, and 39 genes were down-regulated (p-value < 0.05 and |log2(fc)| > 1) (Figures 3C, Supplementary Files 1 and 3). A heatmap and volcano map displayed the distribution of genes in each group ( Figure 3B,D). As demonstrated in Tables 1-3, PrRP stimulated (371 up-regulated genes) more genes than suppressing genes (39 down-regulated genes) on the transcriptional level of BMDMs. In addition, DEGs stimulated by PrRP were composed of the following types of genes: protein_coding genes (94.43%), lincRNA genes (1.03%), antisense genes (0.66%), and other genes (3.88%) ( Figure 3A). Next, the morphological characters of BMDMs before and after PrRP exposure were examined with an optical microscope. As shown in Figure 2C, BMDMs were oval before PrRP (1 nM) exposure, and the shape of BMDMs did not change remarkably after PrRP exposure for 18 h. Overall, PrRP demonstrated no significant effect on the shape of BMDMs.

Identification of DEGs
To explore the effect of PrRP on the transcriptome of BMDMs, cell samples were detected by RNA-seq sequencing ( Figure 3A-D, Table S1, and Supplementary File 2). The quality control test results showed that our RNA-seq was qualified (Supplementary File S2 Figure S1A,B). A total of 410 DEGs were obtained, of which 371 genes were upregulated, and 39 genes were down-regulated (p-value < 0.05 and |log2(fc)| > 1) ( Figure 3C, Supplementary Files 1 and 3). A heatmap and volcano map displayed the distribution of genes in each group ( Figure 3B,D). As demonstrated in Tables 1-3, PrRP stimulated (371 up-regulated genes) more genes than suppressing genes (39 down-regulated genes) on the transcriptional level of BMDMs. In addition, DEGs stimulated by PrRP were composed of the following types of genes: protein_coding genes (94.43%), lincRNA genes (1.03%), antisense genes (0.66%), and other genes (3.88%) ( Figure 3A).

Functional and Pathway Enrichment Analysis of DEGs
To explore the biological meanings of DEGs, several tools were used to analyze the functions of DEGs, such as Metascape, KEGG, and PANTHER.
The Metascape online website was employed to study the functional enrichment of DEGs. All DEGs were subjected to Metascape to acquire the enrichment pathways. The down-regulated DEGs activated several pathways such as the glycosaminoglycan biosynthetic process, brain development, sodium ion transport, and central nervous system neuron differentiation ( Figure 4A,D,G). Meanwhile, enriched pathways stimulated by up-regulated DEGs mainly included the regulation of defense response, response to interferon-beta, response to virus, response to interferon-gamma, inflammatory response, Herpes simplex infection ( Figure 4B,E,H). Besides, all DEGs activated a series of pathways, including regulation of defense response, response to interferon-beta, response to virus, response to interferon-gamma, inflammatory response, and Herpes simplex infection ( Figure 4C,F,I). Finally, these pathways were clustered, followed by connecting with various network diagrams ( Figure 4G-I).
The online tool PANTHER was also employed to investigate the enriched processes of DEGs. All DEGs were classified into the following categories: molecular function (MF), biological process (BP), and cellular compartment (CC).

Protein-Protein Interaction (PPI) Network Analysis
To explore the DEGs-induced protein-protein interactions, DEGs were analyzed with the online tool STRING database, followed by visualizing with Cytoscape. By using the Cytoscape's plug-in cyto-Hubba, a total of eight hub genes (IFIT1, OASL2, IRF7, IFIT3, IFIT2, USP18, IFI44, and RTP4) were acquired with the highest scores. Besides, the Cytoscape plug-in ClueGO was used to study DEGs-induced functional processes.
In addition, the GO analysis and detailed information of eight hub genes were listed in Tables 4 and 5. Besides, the PPI interaction between DEGs was analyzed with the Cytoscape Plug-in Mcode. As demonstrated in Figure 7, four networks were activated by the up-regulated DEGs while no network was stimulated by the down-regulated DEGs.   In addition, the GO analysis and detailed information of eight hub genes were listed in Tables 4 and 5. Besides, the PPI interaction between DEGs was analyzed with the Cytoscape Plug-in Mcode. As demonstrated in Figure 7, four networks were activated by the up-regulated DEGs while no network was stimulated by the down-regulated DEGs.

Verification of Hub Genes with qPCR and Western Blot
In order to verify the accuracy of the RNA-seq data, qPCR was carried out. These hub genes were all protein-coding genes (IFIT1, OASL2, IRF7, IFIT3, IFIT2, USP18, IFI44, and RTP4). As displayed in Figure 8A, PrRP (1 nM) up-regulated the mRNAs of eight hub genes and no hub genes were down-regulated, which were consistent with the RNA-seq results. Moreover, the expression of the hub protein was also detected by the Western blot.
Since only a few commercial antibodies were available for us, the protein expression data of two hub proteins (IFIT1 and USP18) were acquired. As demonstrated in Figure 8B, PrRP showed no noticeable effect on the protein expression of IFIT1 and USP18.

Verification of Hub Genes with qPCR and Western Blot
In order to verify the accuracy of the RNA-seq data, qPCR was carried out. These hub genes were all protein-coding genes (IFIT1, OASL2, IRF7, IFIT3, IFIT2, USP18, IFI44, and RTP4). As displayed in Figure 8A, PrRP (1 nM) up-regulated the mRNAs of eight hub genes and no hub genes were downregulated, which were consistent with the RNA-seq results. Moreover, the expression of the hub protein was also detected by the Western blot. Since only a few commercial antibodies were available for us, the protein expression data of two hub proteins (IFIT1 and USP18) were acquired. As demonstrated in Figure 8B, PrRP showed no noticeable effect on the protein expression of IFIT1 and USP18.
(A) . Verification of hub genes. BMDMs were incubated with PrRP (1nM) for 18 h, followed by a qPCR detection and a Western blot. (A) Total RNA was extracted, and a qPCR detection was carried out to identify eight hub genes. The mRNA level was normalized by the expression of GAPDH. *, significantly different from the control group, ** p < 0.01; *** p < 0.001. Each detection was conducted three times in duplicate. The data were presented as the means ± S.E.M. Statistical significance analysis was performed with the t-test method. (B) BMDMs were lysed and subsequently subjected to Western blot examination (n = 3). The data are shown as the means ± S.E.M. *, significantly different from the control group; ns, no significance. Statistical significance analysis was conducted with the t-test approach.

Protein Modeling of Hub Proteins
To investigate the protein structure of hub proteins, the three-dimensional structures of hub proteins were built using the Modeller (9v23). A total of 1000 models for hub proteins were acquired (except for IFI44 and RTP4), and the model with the lowest DOPE value was collected (IFIT1:  Table 7).

Protein Modeling of Hub Proteins
To investigate the protein structure of hub proteins, the three-dimensional structures of hub proteins were built using the Modeller (9v23). A total of 1000 models for hub proteins were acquired (except for IFI44 and RTP4), and the model with the lowest  Table 7).

Molecular Dynamics Simulation of Hub Proteins
Molecular dynamics (MD) simulation (at least 300 ns) was performed to study the behavior of the hub proteins (Figures 9 and 11). The structural convergence data of hub proteins were displayed as RMSD (Figure 12 Table S3).
To investigate the flexibility of hub proteins, the RMSF of the atoms of hub proteins was analyzed ( Figure 13). Residues with high RMSF suggested high flexibility, whereas low RMSF values indicated few fluctuations between average positions and residues. As displayed in Table S3 and Figure 13

Peptide-Hub Protein Docking
To explore the possible interaction between PrRP and hub proteins, the Rosetta program was employed to analyze the possible dock binding sites for PrRP-hub proteins. As displayed in Figure 15, PrRP was not wholly embedded in the protein structure. As expected, PrRP bound to the N-or C-terminal region of the hub protein, where located at the outside of the protein structure. The radius of gyration (Rg) suggested the compression of the protein atoms, and a decrease in Rg indicated an unstable protein structure. As displayed in Figure 14

Peptide-Hub Protein Docking
To explore the possible interaction between PrRP and hub proteins, the Rosetta program was employed to analyze the possible dock binding sites for PrRP-hub proteins.
As displayed in Figure 15, PrRP was not wholly embedded in the protein structure. As expected, PrRP bound to the N-or C-terminal region of the hub protein, where located at the outside of the protein structure.

Expression of GPR10 on BMDMs
The expression of GPR10 was examined in BMDMs by immunofluorescence stain and Western blot ( Figure 16). As shown in Figure 16C, GPR10 protein was distributed on the cell membrane of BMDMs, and no signal was acquired with the IgG control (the negative control). PrRP 1 nM exposure (18 h) provoked a remarkable increase in the GPR10 protein ( Figure 16A,B) than the control group.

Expression of GPR10 on BMDMs
The expression of GPR10 was examined in BMDMs by immunofluorescence stain and Western blot ( Figure 16). As shown in Figure 16C, GPR10 protein was distributed on the cell membrane of BMDMs, and no signal was acquired with the IgG control (the negative control). PrRP 1 nM exposure (18 h) provoked a remarkable increase in the GPR10 protein ( Figure 16A,B) than the control group. . The data were displayed as the means ± S.E.M. *, significantly different from the control group; * p < 0.05; ** p < 0.01; *** p < 0.001; **** p < 0.0001. Statistical significance analysis was conducted using the one-way ANOVA followed by the Tukey post hoc tests approach. (C) The expression of GPR10 in BMDMs was examined by immunofluorescence staining. The nucleus stained with DAPI was blue, and GPR10 stained with anti-GPR10 antibody was green. Scale bar, 75 μm.

PrRP Modulated Different Functional Enrichment Pathways of BMDMs
In 2012, the Romeroa team discoveries that PrRP promotes the production of proinflammatory cytokines (IL-8, IL-12, IL-1β, and IL-6) in Salmon salar, indicating that PrRP is a local transmitter of the innate immune pathway in leukocytes [13]. In the same line, in this study, PrRP significantly promoted inflammatory immune response pathways, in- . The data were displayed as the means ± S.E.M. *, significantly different from the control group; * p < 0.05; ** p < 0.01; *** p < 0.001; **** p < 0.0001. Statistical significance analysis was conducted using the one-way ANOVA followed by the Tukey post hoc tests approach. (C) The expression of GPR10 in BMDMs was examined by immunofluorescence staining. The nucleus stained with DAPI was blue, and GPR10 stained with anti-GPR10 antibody was green. Scale bar, 75 µm.

PrRP Modulated Different Functional Enrichment Pathways of BMDMs
In 2012, the Romeroa team discoveries that PrRP promotes the production of proinflammatory cytokines (IL-8, IL-12, IL-1β, and IL-6) in Salmon salar, indicating that PrRP is a local transmitter of the innate immune pathway in leukocytes [13]. In the same line, in this study, PrRP significantly promoted inflammatory immune response pathways, including regulation of defense response, response to interferon-beta, response to interferon-gamma, inflammatory response ( Figure 4B,E,H). Overall, these data suggested that PrRP may be involved in inflammation and immune regulation processes.
A series of reports show that PrRP plays a vital role in regulating food intake and energy metabolism [8,[26][27][28][29][30]. In the present study, PrRP regulated several functional pathways responsible for macrophages, including neuron differentiation of the central nervous system, sodium ion transport, brain development, and peptide hormone secretion (Figure 4), suggesting that the expression profile of related genes was affected by PrRP. Given that macrophages are deeply involved in regulating multiple physiological processes of neuroendocrine [31], PrRP may affect the physiological activities of the neuroendocrine field by regulating macrophages.

Common Transcription Factors Tied to PrRP-Regulated DEGs in BMDMs
In order to capture the effects of DEGs on the transcription profile of macrophages at the transcriptome level, DEGs were subjected to TRRUST (version 2) for transcription factor analysis. As demonstrated in Table 6, PrRP stimulated a series of general transcription factors, including Irf1 (interferon regulatory factor 1), Irf8 (interferon regulatory factor 8), and Irf4 (interferon regulatory factor 4). These data implied that PrRP might stimulate the related genes governed by these transcription factors. Besides, PrRP activated several differentiation-related transcription factors, including Stat1 (signal transducer and activator of transcription 1), Nfkb1 (nuclear factor of kappa light polypeptide gene enhancer in B cells 1, p105), Jun (jun proto-oncogene), Ikbkb (inhibitor of kappaB kinase beta), and Stat3 (signal transducer and activator of transcription 3), indicating that PrRP might be involved in the biological processes mastered by above transcription factors. Additionally, the number of down-regulated DEGs is too small to obtain transcription factors.

The Concentration of PrRP in the Experimental System for High-Throughput Sequencing
In this study, the treatment concentration of PrRP on cells is a question worth of concern. To the best of our knowledge, the use of high-throughput sequencing (such as RNA-seq) to investigate the effects of PrRP on the transcriptome of target tissues or cells has not been reported yet. In our experimental system, PrRP 1 nM was selected as the concentration to treat cells based on the following considerations.
On the one hand, the concentration of PrRP was limited to a range as close as possible to the physiological concentration in the body (around nM range), which would be helpful to simulate the real effect of neuropeptides in physiological conditions as much as possible. On the other hand, NPFF, another neuropeptide belonging to the RFamide peptide family as PrRP, provides us with a reference. Recently, Waqas and colleagues investigate the effects of NPFF on the gene expression of J774A.1 macrophages and 3T3-L1 preadipocytes by high throughput sequencing where cells were exposed to NPFF 1 nM for 18 h [15,32]. Moreover, our lab also examine the influence of NPFF on the mouse macrophages RAW 264.7 transcriptome where the same exposure (1 nM, 18 h) was applied to the cells [14]. Therefore, in this study, BMDMs were exposed to 1 nM PrRP for 18 h for subsequent RNA-seq detection.

Expression of PrRP-GPR10 on BMDMs
PrRP-GPR10 is widely expressed in the neuroendocrine system. In the peripheral tissues, as a ligand in the PrRP-GPR10 system, PrRP mRNA is found in the lung, adrenal gland, liver, kidney, pancreas, gut, and reproductive organs [33,34]. In addition, PrRP is also distributed in the NTS, ventrolateral reticular, and DMN nucleus [6]. Meanwhile, as the receptor in the PrRP-GPR10 system, GPR10 mRNA is expressed in the rat adrenal medulla, epididymis, and testis [35]. Given that the pons and hypothalamus are both key parts of controlling food intake and energy metabolism [36], the above data shows that PrRP-GPR10 should play a key role in neuroendocrine metabolism and may even be a potential target for anti-obesity therapy [29,37,38].
In the central nervous system, GPR10 is expressed in the paraventricular hypothalamic (PVN), thalamic reticular (TRN), dorsomedial hypothalamic (DMN) nuclei, periventricular hypothalamic (PEVN), and the nucleus of the solitary tract (NTS) of the brainstem [39]. These data suggest that GPR10 is involved in energy metabolism and food intake [8,40]. Very recently, the Lenka Maletínská team find that GPR10 gene deletion in C57BL/6J mice causes significant metabolic disturbances, as GPR10 KO mice demonstrate enhanced basal neuronal activity, disturbed lipid homeostasis, and altered insulin sensitivity [37].
Moreover, the expression of PrRP-GPR10 in the immune system has already been reported. The Romeroa group demonstrates that PrRP mRNA is expressed in leukocytes from the head kidney and blood of Salmon salar. Moreover, synthetic PrRP provokes the production of pro-inflammatory cytokines (IL-8, IL-12, IL-1β, and IL-6), implying that PrRP could be a local transmitter of the innate immune pathway in leukocytes [13]. In the present study, GPR10 was found to be expressed on the cell membrane of BMDMs ( Figure 16C), which increased after PrRP (1-1000 nM, 18 h) exposure ( Figure 16A,B). These data indicated that PrRP might be involved in diverse physiological processes controlled by macrophages.

Prolactin and Inflammatory Processes
To the best of our knowledge, there are few studies on the involvement of PrRP in the inflammatory process. However, prolactin, secreted by neurons stimulated by PrRP, is widely involved in inflammation and immune processes.
Prolactin is a peptide hormone that is generated in the anterior pituitary gland and in various sites outside of the pituitary. It has been reported that prolactin is involved in a variety of biological processes, including lactation, reproduction, and immune functions. Elevated serum prolactin concentration is often associated with inflammation and immune response [41].
The Jörg-Matthias Brand team finds that prolactin stimulates a pro-inflammatory immune response on peripheral inflammatory cells [42]. Prolactin, at concentrations achievable during medication, pregnancy, and anesthesia, significantly enhances the synthesis of tumor necrosis factor-alpha (TNF-alpha) and interleukin (IL)-12 in lipopolysaccharideactivated human whole blood cultures. These data suggest that prolactin may affect pathophysiological processes in physiological hyperprolactinemic disorders [42].
Actually, prolactin has been considered an inflammatory cytokine, which inhibits the negative selection of autoreactive B lymphocytes [41]. Up to now, hyperprolactinemia has been found in patients with a variety of autoimmune diseases, including systemic lupus erythematosus, rheumatoid arthritis, multiple sclerosis, systemic sclerosis, and autoimmune thyroid disease [41,[43][44][45]. These data indicate that prolactin is involved in the pathological process of the above diseases.

Prolactin and Macrophages
Prolactin and its receptors are widely involved in the functional regulation of macrophages. On the one hand, prolactin promotes the release of cytokines, chemokines, and reactive oxygen species in macrophages [41,[46][47][48][49][50]. On the other hand, the prolactin receptor is widely distributed through the immune system, such as macrophages, lymphocytes, monocytes, granulocytes, natural killer cells, and thymic epithelial cells [51]. Hence, prolactin and its receptor are widely involved in the immune response process.
In the present study, PrRP differentially regulated the gene expression of mouse BMDMs. Besides, GPR10, the receptor of PrRP, was expressed on mouse BMDMs. Given that BMDMs are deeply involved in regulating the neuroendocrine system, our data suggested that PrRP may be involved in various physiological processes controlled by macrophages.

PrRP and Microglia
Recently, the potential neuroprotective effect of PrRP has also been reported. By using a model of AD-like β-amyloid (Aβ) pathology (double transgenic APP/PS1 mice), Lenka Maletínskáa's team find that a PrRP analog (palm 11 -PrRP31) reduces the total amount of senile Aβ plaques in APP/PS1 mice. Moreover, palm 11 -PrRP31 reduces the marker (ionized calcium-binding adaptor molecule 1 (Iba1)) of microglia near the Aβ plaque, suggesting that PrRP may have a potential neuroprotective effect [52].
Microglia are a type of cells derived from bone marrow hematopoietic stem cells. In the early stage of embryonic development, bone marrow hematopoietic stem cells enter the central nervous system and finally differentiate into microglia through the monocyte lineage [53]. As the resident macrophages in the central nervous system (CNS), microglia play an essential role in the innate and acquired immune response in local regions [54]. Given that microglia are mainly distributed in large non-overlapping areas throughout the CNS [55,56], and these areas may overlap with areas where PrRP is distributed (such as pons and hypothalamus) [33][34][35]40]. Therefore, a question arises: does PrRP regulate glial cells? Unfortunately, to the best of our knowledge, the research on the effect of PrRP on microglia has not been reported yet, which will be an exciting direction worth exploring.

Considerations of Double Positive Cells in the Control Group of Flow Cytometry Results
In the present study, mouse BMDMs were differentiated from primary mouse bone marrow cells induced by cytokines. It is worth noting that the results of the control group (no antibody was added) showed that 0.14% of the cells were still positive for anti-F4/80 and anti-CD11b ( Figure S1A). Regarding this phenomenon, we believe that the following reasons may be responsible for this.

Non-Specific Fluorescent Signal Caused by Dead Cells
Primary cells are of great significance in biomedical research because they have certain advantages over cell lines in maintaining the physiological environment of biological samples. However, a significant shortcoming in primary cell research is the difficulty in obtaining a single type of cell with 100% purity. In our experimental system, the monocytes in the bone marrow of mice were isolated according to a widely used method [57][58][59], followed by cytokine induction and eventually formed BMDMs. It should be pointed out that this method has a certain probability of introducing non-specific cells that initially existed in the bone marrow of mice. These cells may die in an unsuitable environment and eventually introduce a non-specific fluorescent signal in the flow cytometry test results. Moreover, the fluorescent signal of dead cells is often displayed as a diagonal signal on the scatter plot of a flow cytometer (it should be noted that the fluorescent signal on the diagonal of the scatter plot does not mean that these signals are certainly come from dead cells, as these signals may be a superposition of signals from live and dead cells) [60][61][62]. Furthermore, the intensity of non-specific fluorescence produced by dead cells can be similar to the intensity of fluorescence produced by fluorescein, which may introduce non-specific fluorescence signals into the final flow cytometer results [60][61][62]. In this study, the fluorescence signal of the control group (no antibody was added) was on the diagonal line ( Figure 2B). Therefore, there is a possibility that the double-positive signal (only 0.14%) in the control cells may be part of the fluorescent signal introduced by dead cells.

Non-Specific Fluorescent Signal from Cell Debris or Tiny Tissue Pieces
Cells or cell-like particulate matter is the object of flow cytometry analysis, and cell debris is inevitably present in the process of flow cytometry detection [62]. In this study, a series of measures were adopted to reduce non-specific signals detected by flow cytometry.
Firstly, a widely used sterile filter (100 mesh) was used to obtain a single-cell suspension as much as possible during the experiment of obtaining mouse primary bone marrow cells. Secondly, in the flow cytometry detection step, a routinely used threshold was set to eliminate obvious cell debris and other non-target cells as much as possible.
Even so, our experimental system was still unable to absolutely distinguish the tiny amounts of other components that may be introduced during the isolation and culture of primary cells from BMDMs (especially those with tiny tissue masses similar in size to BMDMs). These non-specific substances might cause minor non-specific signals in the final results. In our experimental system, the double-positive rate of BMDMs in the experimental group was 96.7%, which was significantly higher than the control group (no antibody was added, double-positive rate: 0.14%), indicating that the purity of BMDMs was acceptable.

The Excitation Wavelength of Fluorescein Is the Same as the Emission Wavelength, Which Is a Wavelength Range, Rather Than a Specific Value
In this study, PE-conjugated anti-F4/80 and FITC-conjugated anti-CD11b were used. In addition, the detection results of the flow cytometer have also undergone channel compensation processing as usual [63][64][65]. However, the excitation wavelength of fluorescein (the same as the emission wavelength) belongs to a wavelength range rather than a specific value [64]. Hence, it is difficult to obtain a completely 100% pure signal. In our opinion, the above factors may also contribute to the small amount of double-positive fluorescence signal in the control group of the flow cytometry test result ( Figure 2B).

Non-Specific Fluorescence Signals Caused by BMDMs
The intensity of the non-specific fluorescence of the cell is determined by the cell itself and is related to the size of the cell. Generally speaking, non-specific fluorescence produced by large cells is strong, while non-specific fluorescence produced by small cells is weak. Larger macrophages have stronger non-specific fluorescence than smaller lymphocytes [62,63,66]. In our study, flow cytometry was used to detect the purity of BMDMs, and the large volume of BMDMs might also introduce non-specific fluorescent signals for the control group ( Figure 2B).

The Limitations of Our Study
Firstly, the changes in the protein levels of hub proteins are worth studying. In this study, through a series of bioinformatics methods, eight hub genes were obtained, which may be the critical nodes of PrRP regulating the gene network of BMDMs. Although the RNA-seq results of the hub genes were verified by qPCR, the changes in the protein levels of the eight hub genes would provide a solid biological basis for the significance of this study. However, limited by the fact that only a few commercial antibodies were currently available to us, the two hub proteins' Western blot data (IFIT1 and USP18) were obtained. As shown in Figure 8B, PrRP did not significantly change the protein levels of IFIT1 and USP18 ( Figure 8B). In our opinion, the following concerns would be helpful to understand this phenomenon.
(1) Hub genes responded to PrRP stimulation in protein levels might be late than in mRNA level. After PrRP (1 nM) treatment for 18 h, the mRNA levels of IFIT1 and USP18 were significantly increased. However, the influence of PrRP on the protein expression of hub genes might be longer than 18 h. In the process of gene expression (from mRNA to protein), it takes minutes to hours to translate mRNA into proteins [67]. (2) PrRP may regulate the expression of the hub proteins of BMDMs in a variety of ways. Although PrRP (1 nM) treatment for 18 h significantly affected the expression of hub genes, this does not mean that PrRP will cause changes in hub proteins' expression levels. In the step from mRNA to protein, a variety of post-transcriptional modifications can cause changes in protein levels, including mRNA degradation control, mRNA translocation control, protein degradation control, translation control, mRNA translocation control, and protein degradation control [67]. Therefore, these above links may be employed by PrRP to regulate the expression of the hub proteins of BMDMs. (3) PrRP may modulate the functions of the hub proteins of BMDMs in diverse manners.
It is the function of proteins, rather than the protein expression, that plays a pivotal role in cellular activities [68]. Various biologically active molecules (including neuropeptides) may regulate each hub protein's activity in a variety of ways, including phosphorylation, heterogeneous regulation, covalent modification, etc. [67]. However, this study's focus is to investigate the effect of PrRP on the transcriptomic gene expression of BMDMs, which may provide clues for the subsequent exploration of the functional regulation of PrRP on BMDMs. Therefore, follow-up studies on the complex biological functions of PrRP on BMDMs are worth looking in to. (4) The regulation of PrRP on the transcriptome gene expression of BMDMs may be more complicated than we previously assumed. Since proteins and various proteinprotein interaction (PPI) networks are prominent members that play an essential role in regulating various biological processes in cells, the influence of neuropeptide PrRP on BMDMs may be a complicated process. In the same vein, emerging data have indicated that PrRP may have a wide range of effects in regulating the neuroendocrine system [40,69].
Secondly, the binding mode of PrRP and hub proteins needs to be verified by experiments. In this study, with the help of molecular simulation tools, the possible binding modes of PrRP and hub proteins were predicted, which may provide clues for exploring the mode of action of PrRP. However, reliable experimental verification is necessary to understand the mechanism of PrRP fully.
Finally, the key "driver" of the network formed by PrRP-activated DEGs needs to be further studied. In the present study, our data described the gene-expression signature of murine bone marrow-derived macrophages from the transcriptome level, which may provide preliminary clues for the understanding of the possible effects of PrRP on the immune system. However, considering the complexity of the gene expression regulatory network [70], it is of great significance to identify the essential molecules of PrRP to regulate the gene expression of BMDMs from the experimental level-such as these "driver" genes.

Ethical Statement
The study was conducted according to the guidelines of the Declaration of Helsinki, and approved by the Ethics Committee of Northwestern Polytechnical University (protocol code 201900048, 2 January 2020).

Animals
Male mice (C57BL/6 strain, 18-22 g) were kept in approved plastic cages with humidity of 65-74% and temperature of 20-22 • C and a photocycle of 12 h dark/light. The mice had free access to water and food. All measures were conducted to minimize animals suffering during the whole experiment process.
Cell culture dishes were purchased from Corning, Inc. (Corning, NY, USA). Prime-Script TM 1st Strand cDNA Synthesis Kit and SYBR ® Premix Ex Taq TM II kit were acquired from TaKaRa (Dalian, China). QiaQuick PCR extraction kit was obtained from Qiagen (Venlo, The Netherlands). The red blood cell lysing solution was from Beyotime (Shanghai, China). L-929 cells were purchased from the Stem Cell Bank of the Chinese Academy of Sciences (Shanghai, China).
BCA Protein Assay Kit was from Thermo Scientific Pierce (Thermo Fisher scientific, MA, USA). PVDF membrane, ECL detection kit, and protease inhibitor cocktail III (EDTAfree) were provided by Millipore Corporation (Bedford, MA, USA). The Rabbit anti-GPR10 polyclonal antibody and horseradish peroxidase conjugated-goat anti-rabbit IgG was acquired from Thermo Fisher Scientific (Beverly, MA, USA). Rabbit anti-IFIT1 monoclonal antibody, rabbit anti-USP18 monoclonal antibody, rabbit anti-Actin monoclonal antibody, and anti-rabbit IgG (Alexa Fluor 488 Conjugate) were obtained from Cell Signaling Technology (Beverly, MA, USA). FITC Rat anti-mouse CD11b antibody and PE Rat-anti mouse F4/80 antibody were from BD Pharmingen (San Diego, CA, USA).
PrRP was synthesized by the GL Biochem Ltd. (Shanghai, China) using the solidphase peptide synthesis method. The mass of the peptide was confirmed using a mass spectrometer (LCMS-2010EV, Shimadzu, Japan). PrRP31 was purified by HPLC, and peptide (purity > 98%) was used in this study.

BMDM Preparation
BMDMs were isolated by flushing mouse tibias and femurs with cold PBS. Cells were centrifuged at 1200 rpm/min for 7 min (4 • C), and the pellet was resuspended with 3 mL red blood cell lysing solution for 5 min. Then, the cell suspensions were centrifuged at 800 rpm/min for 5 min (at 4 • C), and the pellet was resuspended with 10 mL PBS. Next, the cell suspensions were filtered through a sterile filter (100 mesh) to collect a single cell suspension.
The single-cell suspensions were maintained in DMEM complete culture medium (streptomycin (100 µg/mL)/penicillin (100 units/mL), 10% FBS) at 37 • C with 5% CO 2 in a humidified incubator. Fourteen h later, the non-adherent cell supernatant was obtained, centrifuged at 800 rpm/min for 5 min (4 • C), and the pellet was resuspended in DMEM complete cell culture medium. Next, cells were maintained in a 100-mm culture dish supplemented with 20% L-929 conditioned media for 8 d (the fresh medium was refreshed every 3 d). At 8 days, the cells were differentiated into BMDMs, where more than 96% of the cells were double positive for F4/80 and CD11b.

RNA-seq Sample Preparation and Microscope Detection
BMDMs were exposed to PrRP (1 nM) for 18 h, followed by RNA sequencing examination. Cells were washed with PBS (3 times) and lysed with TRIzol (1 mL). Next, cell lysates were quickly stored in liquid nitrogen. Subsequently, the RNA-seq detection was performed by the Novogene Co Ltd. (Beijing, China).
Cell images were acquired by a confocal microscope (Leica TCS SP5, Leica Microsystems, Mannheim, Germany). In addition, the detailed structure of BMDMs was detected by a transmission electron microscope HT7700 (Hitachi High-Technologies, Tokyo, Japan).

Flow Cytometry Detection
Cells were rinsed twice with pre-cooled buffer (BSA-PBS-1%) and washed with trypsin (0.25%). The cells were resuspended in a fresh medium to form a uniform cell suspension, followed by centrifugation (1000/rpm, 5 min). Next, the pelleted cells acquired with centrifugation were resuspended in a buffer solution (BSA-PBS-1%) to form a uniform cell suspension. Five min later, the cells were centrifuged (1000/rpm, 5 min), and the pellet was resuspended in a buffer solution (BSA-PBS-1%). Then, the number of cells were counted (1 × 10 6 /100 µL). The antibody was added into the cell suspension, followed by incubating in the dark for 25 min. During the whole incubation process, the cell mixture was mixed every 2 min to make sure the antibody and the cells to be thoroughly combined. Subsequently, the cell mixture was subjected to centrifugation (1000 rpm/min, 5 min), and the supernatant was discarded. Then, 400 uL of buffer solution (BSA-PBS-1%) was added into the tube to resuspend the cell pellet. Next, the whole-cell suspension was filtered with a 300-mesh sterile filter to ensure that only single cells were subjected to flow cytometry detection. Finally, the purity of macrophages was detected by flow cytometry (BD Calibur, Biosciences, San Jose, CA, USA), followed by analysis with the Cellquest (BD) and Modfit software. Low-quality reads were removed, and the clean reads with high quality were used. The paired-end clean reads were aligned to the musculus reference genome with the Hisat2 v2.0.5. Finally, the FPKM (fragments per kilobase of exon model per million reads mapped) of each gene was calculated, and the read was mapped to each gene by using the FeatureCounts v1.5.0-p3.

Differential Expression Interpretation
Differential expression gene (DEGs) detection was conducted with the DESeq2 R package (1.16.1). The p-values were adjusted with the Benjamini and Hochberg's method to analyze the false discovery rate (FDR). Genes with an adjusted p-value < 0.05 were considered as DEGs. The heat map of DEGs was generated by the toolkit TBtools [71].

GO and KEGG Enrichment Interpretation
In order to analyze the ontology (GO) enrichment of differentially expressed genes (DEGs), the R package "ClusterProfiler" [72] was utilized. GO terms (adjusted p < 0.05) were considered as significantly enriched.
To explore the pathways associated with DEGs, the enrichment study was performed with Metascape online tool (http://metascape.org/, date of retrieving data: 8 November 2020) [73,74], and the significant biological processes were obtained.
To further investigate the functions of DEGs, DEGs were subjected to further analysis using PANTHER (http://www.pantherdb.org/, date of retrieving data: 11 November 2020) [75]. The main GO categories were acquired, including molecular function (MF), biological process (BP), and cellular component (CC).
Besides, KEGG (http://www.genome.jp/kegg/, date of retrieving data: 11 November 2020) was used to explore the functional enrichment of DEGs. To show the KEGG pathway maps clearly, a Cytoscape plug-in KEGGParser was employed to show biological networks.

Protein-Protein Interaction (PPI) Network Analysis
In order to explore the interaction among DEGs, all DEGs were analyzed with the online tool STRING (http://string-db.org, date of retrieving data: 15 November 2020) [76]. Only experimentally proved interactions (a combined score > 0.4) were considered significant. Next, the PPI network was built with the Cytoscape software (Ver3.8.0), and the Cytoscape plug-in (Molecular Complex Detection, MCODE) was utilized to analyze the modules of the PPI network. Subsequently, function enrichment analysis for DEGs of the modules was conducted (p < 0.05 was considered as statistically significant).
In order to identify the top-ranked hub genes from the PPI network, CytoHubba (a Cytoscape plug-in) was used in the following analysis. A total of eight hub genes were obtained according to the scores of the following methods, including MCC, MNC, Radiality, DMNC, Degree, Betweenness, Closeness, BottleNeck, EcCentricity, EPC, Stress and Clustering Coefficient [79,80].

qPCR Analysis
This approach has been described elsewhere [81]. Briefly, total RNA was extracted from cells with the TRIzol following the manufacturer's instruction. Then, cDNA was transcribed from 1 µg of RNA with a PrimeScript TM 1st Strand cDNA Synthesis Kit. Next, the gene expression level was detected with the SYBR ® Premix Ex Taq TM II system and the MX3000P Real-Time PCR System (Stratagene). Real-time PCR conditions were set as follows with minor modifications: 94 • C for 30 s, 95 • C for 5 s, 56 • C. Data were normalized by GAPDH data using the comparative 2 −∆∆CT approach. Primers were shown in Table 9. Data were acquired in duplicates three times.

Western Blot Detection
This approach has been previously described [82]. In brief, cells were washed with phosphate-buffered saline3 times and lysed by lysis solution (150 mM NaCl, 5 mM EGTA, 1% Nonidet P-40, 50 mM Tris/HCl, 0.5% sodium deoxycholate, 1 unit protease inhibitor cocktail III (EDTA-free), 0.1% SDS, PH 7.4). Then, cell lysates were centrifuged (12,000× g for 12 min, 4 • C), and the protein concentration of the supernatants was evaluated by the BCA Protein Assay Kit. Samples for immunoblot (18-24 µg/lane) were analyzed with the Bio-Rad mini-gel system by 10% SDS-polyacrylamide gel electrophoresis. Next, the proteins were blotted onto PVDF membranes using the Bio-Rad wet blotter system. The membranes were blocked with the solution (5% non-fat milk in the tris-buffered saline containing 0.05% tween-20) for 1 h. Then, the membranes were rinsed three times with TBST and were treated at 4 • C with appropriate antibodies (anti-GPR10 was 1:1000, and anti-Rabbit antibody was 1:10,000, and anti-Actin was 1:5000) for 12 h. The membranes were rinsed with TBST 3 times, followed by treated with horseradish peroxidase-conjugated secondary antibodies at room temperature for 2 h. Subsequently, the PVDF membranes were detected using an ECL detection kit. The band intensity was analyzed by the ChemDocTM XRS (Bio-Rad, Hercules, CA, USA) and Image J [83].

Immunofluorescence Stain
Immunofluorescence stain was conducted as previously described [84]. In brief, BMDMs were cultured on the glass slides, followed by fixing with 4% paraformaldehyde for 15 min. Then, cells were rinsed with PBS 3 times, exposed to 0.1% Triton X-100 (12 min), and treated with 5% normal rabbit serum (3 h). Next, BMDMs were incubated with rabbit anti-GPR10 polyclonal antibody (1:250) at 4 • C for 14 h, followed by incubation with anti-rabbit IgG (Alexa Fluor 488 Conjugate) for 1.5 h. Subsequently, cells were treated with DAPI for 5 min to stain the nucleus. All figures were captured with a confocal microscope (Leica TCS SP5, Leica Microsystems, Mannheim, Germany).

Homology Modeling of Proteins
The three-dimensional (3-D) structures of proteins were built using homology modeling (the modeling of Ifi44 and Rtp4 was conducted with Robetta de novo structure prediction program [85,86] due to the lack of sufficient templates). Briefly, protein sequences were acquired from the NCBI nucleotide database (https://www.ncbi.nlm.nih.gov/protein/, date of retrieving data: 27 September 2020), and the blast module was employed to align the hub protein sequence in the PDB database [87]. Protein templates were retrieved from the online RCSB PDB Protein Data Bank. Firstly, three templates (query cover > 45%) for each protein were selected for homology modeling. Next, 3-D models of the hub protein were constructed using Modeller (9v23) [88]. Modules for multiple template modeling were used, including Align2d, Salign, and Model. One thousand candidate three-dimensional models for each protein were built, and the model with the lowest discrete optimized protein energy (DOPE) value was finally selected. Finally, all models were subjected to the online tool MolProbity (http://molprobity.biochem.duke.edu/, date of retrieving data: 13 December 2020) for quality evaluation [89].

Molecular Dynamics (MD) Simulation
The molecular dynamics simulation of the hub protein was performed by using the Gromacs 2018.12 [90]. All simulations were carried out in the CHARMM36 force field [91].

Molecular Dynamic Simulation: Protein in Water
Molecular dynamics simulation was performed in a condition with minor modifications. The protein model was solvated in an octahedron box with a TIP3P water model (1.0 nm). The simulated system was neutralized by adding Cl − or Na + ions, and periodic boundary conditions (PBC) were used in all directions. Next, energy minimization for the protein model was conducted with the steepest descent (50,000 steps) with the max force (<100 KJ/mol). Subsequently, the whole simulation was performed under equilibration phases with NVT (100 ps, 298.15 K) and NPT (300 ps, 298.15 K, 1.0 bar), respectively. The whole simulation was conducted at least 300 ns for each protein.

Molecular Dynamic Simulation: Analysis
The MD trajectory was analyzed using GROMACS utilities to obtain the RMSF (root mean square fluctuation), RMSD (root mean square deviation), and radius of gyration (Rg). The 3-D structures of proteins were drawn with the Pymol software (Delano, W.L. The Pymol Molecular Graphics System (2002) DeLano Scientific, SanCarlos, CA, USA. http://www.pymol.org, date of retrieving data: 11 November 2020).

Dock
The 3-D model of PrRP was constructed by PEP-FOLD3 [92], and the 3-D structure of hub proteins was obtained from the MD-optimized protein structures.
(1) The primary docking complex, consisting of PrRP (ligand) and hub protein (receptor), was constructed by ZDOCK (3.02) [93][94][95]. (2) The primary docking complexes were sent to the Flexible peptide docking module of the Rosetta program (3.9) for subsequent docking. (A) Pre-pack mode: 1 model was constructed. (B) Low-resolution ab-initio analysis: 100 models were obtained, and 1 model with the best docking score was chosen for subsequent analysis. (C) Refinement analysis: 100 docking models were acquired, and a docking model was finally selected based on the total_score.

Statistical Analysis
Data were demonstrated as means ± standard error of the mean (S.E.M). Data were interpreted by one-ANOVA followed by the post hoc tests (Tukey). The statistical analysis was performed by using the GraphPad Prism software version 8.0 (San Diego, CA, USA).

Conclusions
Our work showed that PrRP significantly changed the transcriptome profile of BMDMs, implying that PrRP may be involved in various physiological activities mastered by macrophages. Our data provided clues for in-depth exploration of the role of PrRP in regulating diverse physiological processes governed by macrophages.