Effects of interleukin 1β on long noncoding RNA and mRNA expression profiles of human synovial fluid derived mesenchymal stem cells

Synovial fluid-derived mesenchymal stem cells (SFMSCs) play important regulatory roles in the physiological balance of the temporomandibular joint. Interleukin (IL)-1β regulates the biological behavior of SFMSCs; however, the effects of IL-1β on long noncoding RNA (lncRNA) and mRNA expression in SFMSCs in the temporomandibular joint are unclear. Here, we evaluated the lncRNA and mRNA expression profiles of IL-1β-stimulated SFMSCs. Using microarrays, we identified 264 lncRNAs (203 upregulated, 61 downregulated) and 258 mRNAs (201 upregulated, 57 downregulated) that were differentially expressed after treatment with IL-1β (fold changes ≥ 2, P < 0.05). Kyoto Encyclopedia of Genes and Genomes pathway analysis found that one of the most significantly enriched pathways was the NF-κB pathway. Five paired antisense lncRNAs and mRNAs, eight paired enhancer lncRNAs and mRNAs, and nine paired long intergenic noncoding RNAs and mRNAs were predicted to be co-expressed. A network constructed by the top 30 K-score genes was visualized and evaluated. We found a co-expression relationship between RP3-467K16.4 and IL8 and between LOC541472 and IL6, which are related to NF-κB pathway activation. Overall, our results provide important insights into changes in lncRNA and mRNA expression in IL-1β-stimulated SFMSCs, which can facilitate the identification of potential therapeutic targets.

www.nature.com/scientificreports/ LncRNA and mRNA co-expression network construction. The "WGCNA" package in RStudio was used to analyze the lncRNA and mRNA co-expression network 29 . The Soft Threshold value was set as 18, and the "one-step network construction and module detection" method was applied. Differentially expressed lncRNAs and mRNAs showing a fold change > 1 were analyzed. K-score, which represents the central location of a gene within a network, was used to determine the crucial genes in the co-expression network.

Reverse transcription quantitative PCR (RT-qPCR).
Five randomly differentially expressed lncRNAs and mRNAs were analyzed by qPCR. We obtained the cDNA using a Superscript Reverse Transcription system (Takara Bio, Shiga, Japan). Relative gene expression was determined by qPCR using the SYBR Green qPCR Kit (Roche, Basel, Switzerland) and the Roche LightCycler 480 System (Roche). The 2 −ΔΔCt method was used to calculate gene expression levels. Primers used are listed in Table S1. The primer sequences were designed using PrimerExpress software (Applied Biosystems, Foster City, CA, USA) by a technologist at Generay Biotech Co., Ltd., Guangzhou, China and was validated in Pubmed. Primers were synthesized by RiboBio Co., Ltd., Guangzhou, China. The unique products after PCR were confirmed by melt-curve analysis using the LightCycler 480.
Western blot analysis. Cells were lysed in RIPA Lysis Buffer (P0013B, Beyotime, China) supplemented with Protease Inhibitor Cocktail (CW2200, CWBIO, China) and phosphatase inhibitor (CW2383S, CWBIO, China) to isolate proteins, whose concentrations were determined using BCA assays (CW0014S, CWBIO, China). The proteins were then incubated at 99 °C for 10 min for denaturation and separated using 6-8% sodium dodecyl sulfate-polyacrylamide gel electrophoresis (P0012A, Beyotime, China), followed by transfer to polyvinylidene fluoride membranes (ISEQ00010, Millipore, MA). The membranes were blocked with 5% nonfat milk (232100, BD Biosciences, CA) for 1 h at room temperature. The blots were then cut and incubated with the following primary antibodies: rabbit anti-phospho-p65 (  Immunofluorescence assay. Cells cultured in confocal dishes were fixed in 4% paraformaldehyde for 15 min. Subsequently, the cells were washed with PBS three times and treated with 0.3% Triton for 15 min. The cells were then incubated in 5% nonfat milk for 1 h at room temperature. Anti-NF-κB p65 rabbit monoclonal antibodies (4764S, 1:1000; RRID: AB_10859369; Cell Signaling Technology, MA) were added to the dishes, and the cells were incubated overnight at 4 °C. PBS was added to one dish as a negative control. The secondary antibody DyLight 594 AffiniPure Goat anti-rabbit IgG (E032320-01; 1:100; EarthOx, USA) was then added, and plates were incubated in the dark for 1 h. Plates were washed with PBS three times, and nuclei were counterstained with 4′,6-diamidino-2-phenylindole for 5 min. Finally, the cells were imaged using a laser-scanning confocal microscope (Carl Zeiss, Oberkochen, Germany). Statistical analysis. Differential expression levels of lncRNAs and mRNAs were analyzed by Student's t tests between two groups. Differences between more than two groups were analyzed by ANOVA and Tukey test. Fisher's exact tests were used for GO and KEGG pathway analysis. Statistical analysis was conducted using SPSS 22.0 (IBM, Armonk, NY). Results with P values of less than 0.05 were considered significant.

Results
LncRNA and mRNA expression changes in SFMSCs after IL-1β stimulation. SFMSCs were isolated from synovial fluid samples which were donated by 3 temporomandibular joint osteoarthritis patients when they accepted arthrocentesis and lavage treatment. Differentially expressed LncRNAs and mRNAs in SFMSCs before and after treatment with IL-1β in vitro were detected using human Arraystar Human LncRNA Microarray V3.0. In total, we identified 14,352 lncRNAs varied slightly in expression after treatment with IL-1β in SFMSCs, including 8014 upregulated and 6338 downregulated transcripts. Of the 20,332 mRNAs varied slightly in expression, 8382 were upregulated and 11,950 were downregulated (fold change cut-off = 1.0) (Fig. 1A, B; Table S2). Among the target molecules, 201 mRNAs and 203 lncRNAs were upregulated, whereas 57 mRNAs and 61 lncRNAs were downregulated (fold changes ≥ 2, P < 0.05) (Fig. 1C, D; Table S4). There were 27 upregulated mRNAs, 19 upregulated lncRNAs, and only one downregulated lncRNA showing fold changes ≥ 10 in SFMSCs after IL-1β stimulation (P < 0.05) ( Fig. 1C-F). Only 41 mRNAs and 24 lncRNAs were significantly upregulated, whereas 6 mRNAs and 1 lncRNAs were significantly downregulated among the target molecules (fold changes ≥ 2, P < 0.05, FDR < 0.05) ( Table S5). The distribution of differentially expressed lncRNAs and mRNAs in chromosomes is shown in Fig. 2A. We identified 262 well-annotated lncRNAs with fold changes ≥ 2; these lncRNAs were classified into six categories: 16 were bidirectional, 25 were exon sense-overlapping, 121 were intergenic, 38 were intron sense-overlapping, 29 were intronic antisense, and 33 were natural antisense RNAs (Fig. 2B). We then picked out the antisense lncRNAs (corresponding to the intronic antisense lncRNAs and natural antisense lncRNAs in Fig. 2B), enhancer lncRNAs which were reported by U. A. Orom et al. 30 , and long intergenic noncoding RNAs (lincRNAs, corresponding to the intergenic lncRNAs in Fig. 2B) and their associated mRNAs. Those lncRNAs and their associated mRNAs with fold changes ≥ 2 and P < 0.05 are shown in Fig. 2C-E.

GO enrichment and KEGG pathway analysis of differentially expressed mRNAs in SFMSCs after IL-1β stimulation.
We subjected differentially expressed mRNAs showing fold changes ≥ 2 and P values < 0.05 to GO analysis. Significantly upregulated mRNAs were mainly enriched in GO terms belonging to biological processes and molecular functions, with one GO term belonging to cellular components (IκB/ NF-κB complex; Fig. 3A-C). In contrast, the significantly downregulated mRNAs were only enriched in 3 GO terms belonging to the biological process category and 2 GO terms belonging to the molecular functions category (Fig. S1). According to KEGG pathway analysis, only the significantly upregulated mRNAs were enriched, whereas the downregulated mRNAs showed no pathway enrichment meeting the screening criteria (P < 0.05, q < 0.05). The significantly upregulated mRNAs were mainly enriched in TNF signaling pathway, cytokinecytokine receptor interaction, and NF-κB signaling pathway according to the KEGG pathway analysis (Fig. 3D). We then visualized the upregulated mRNAs enriched in NF-κB signaling pathway (Fig. 3E). We also visualized the network constructed by genes belonging to the top three most significantly enriched pathways (Fig. 3F).
To further investigate mRNA functions at the protein level, we analyzed the interactions of proteins encoded by differentially expressed mRNAs showing fold changes ≥ 2 in the STRING database. The network contained 98 nodes and 209 edges (Fig. S2A); disconnected nodes in the network were hidden. Next, K-scores were used to evaluate the importance of genes; the network constructed by the top 30 K-score genes was visualized (Fig. S2B) and found to consist of two core subnetworks (Fig. 4A, B).
LncRNA/mRNA co-expression network in SFMSCs after IL-1β stimulation. The functions of most lncRNAs have not yet been annotated. Therefore, we analyzed lncRNA functions by evaluating the co-expressed mRNAs. The co-expression profiles of differentially expressed lncRNAs and mRNAs showing fold changes > 1 were analyzed using weighted correlation network analysis. The genes were subdivided into 30 co-expression modules. We found that only MEgreen and MEcyan modules was obviously positive correlated with the treating factor (IL-1β, P value < 0.05, the genes in these 2 modules were listed in Table S6), and the MEgreen module was most significantly related to IL-1β stimulation (r = 0.88, P = 0.02; Fig. 5A). We then identified the upregulated mRNAs enriched in NF-κB signaling pathway (in Fig. 3F) and their co-expressed lncRNAs and mRNAs (weight value ≥ 0.4) in MEgreen module and then visualized the network constructed by these genes. The network contained 131 nodes and 572 edges (Fig. 5B). The most crucial subnetwork was subsequently constructed in Cytoscape and was found to contain 15 nodes and 39 edges. LncRNA AK022421, BC046172, AF052103, BU564349 and AA992250 were predicted to be associated with the NF-κB signaling pathway (Fig. 5C). The network constructed by genes showing the top 1000 weight values in the MEgreeen module was then visualized, and was found to contain 150 nodes and 1000 edges (Fig. 6A). The most crucial subnetwork was constructed by genes with the highest K-scores, including three lncRNAs (RP3-467K16. 4 www.nature.com/scientificreports/ mRNAs (ICAM1, IL-8, and ACP2; Fig. 6B). RP3-467K16.4 showed the highest K-score in the network shown in Fig. 6B. The 39 mRNAs in Fig. 6A showing co-expression with RP3-467K16. 4 were selected, and the network was reconstructed (Fig. 6C). According to KEGG pathway analysis, we found that RP3-467K16.4 may be related to pathways contributing to inflammation activation, such as the tumor necrosis factor and NF-κB pathways (Fig. 6D).
Co-expression relationship between RP3-467K16.4 and IL8 and between LOC541472 and IL6, which related to NF-κB pathway activation. The NF-κB pathway was activated after IL-1β stimulation according to immunofluorescence results, consistent with our KEGG analysis results (Fig. 7A, B). To verify that the lncRNAs predicted in Fig. 2C-E were related to their associated mRNAs, we selected one lncRNA (LOC541472) and associated mRNA (IL-6) pair and analyzed their relationship in SFMSCs after IL-1β stimulation. We found that both IL-6 and LOC541472 were significantly upregulated after IL-1β stimulation. When we knocked down LOC541472, both gene and protein expression levels of IL-6 and the NF-κB pathway were downregulated significantly (Fig. 7C-F, H, I). We detected strong LOC541472 expression in the nucleus, but weak IL-6 fluorescence in the control group. After IL-1β stimulation, both LOC541472 and IL-6 translocated into the cytoplasm. Interestingly, the localization of LOC541472 and IL-6 in SFMSCs was consistent with the RNAscope results (Fig. 7G). To verify the co-expression relationships predicted by the weighted group correlation network analysis, we selected the lncRNA with the highest K-score (RP3-467K16.4) and IL-8 mRNA, which acts downstream of the NF-κB pathway, in the subnetwork shown in Fig. 6C. We found that both RP3-467K16.4 and IL-8 . (E) Differentially expressed enhancer-like lncRNAs and their nearby coding genes (distance < 300 kb); both the lncRNA and its paired mRNA showed fold changes ≥ 2 and P < 0.05 among all differentially expressed lncRNAs and mRNAs. Diamonds represent lncRNAs, and circles represent mRNAs; red represents upregulation, and blue represents downregulation; different sizes represent the absolute value of the fold change. Protein names were used to label the genes in (A); Genesymbol was used to label the genes in (C-E). Images (C-E) were created using CytoScape 3.8.0 Software. exon sense-overlapping: the LncRNA's exon is overlapping a coding transcript exon on the same genomic strand; intronic: the LncRNA is overlapping the intron of a coding transcript on the same genomic strand; natural antisense: the LncRNA is transcribed from the antisense strand and overlapping with a coding transcript; non-overlapping antisense: the LncRNA is transcribed from the antisense strand without sharing overlapping exons; bidirectional: the LncRNA is oriented head to head to a coding transcript within 1000 bp; intergenic: there are no overlapping or bidirectional coding transcripts nearby the LncRNA. www.nature.com/scientificreports/ were significantly upregulated after IL-1β stimulation. Additionally, knockdown of RP3-467K16.4 downregulated expression of the IL-8 gene and components of the NF-κB pathway (Fig. 7J-O).

Discussion
The etiology of temporomandibular joint disorders is complex and unclear. The inflammatory factor IL-1β is known as an important pathogenic factor for temporomandibular joint osteoarthritis 5 . One of the most important underlying pathogenic mechanisms involves the ability of IL-1β to affect joint injury repair by inhibiting the chondrogenic potential of MSCs; NF-κB pathway activation and upregulation of inflammatory cytokines, such as IL6, contribute to this process 4,6,31,32 . In addition, IL-1β can also affect the immune regulation function of MSCs in the joint. During the early phase of inflammation in the synovium, the resident macrophages secrete a large number of inflammatory factors, such as IL-1β, which shifts the synovial tissue from homeostasis toward catabolism, promoting joint destruction. The macrophages also promote CXCR3 + CD4 + Th1 cell migration into the synovium, which promotes polarization of resident macrophages toward M1 phenotype and enhances  www.nature.com/scientificreports/  www.nature.com/scientificreports/ inflammatory factor secretion, creating a "positive inflammatory loop" that results in destruction progression 33,34 . The cartilage damage caused by the inflammatory immune response requires MSCs to spontaneously differentiate into chondrocytes for repair 35 . Meanwhile, the upregulated IL-1β enhances IL-1Ra secretion of MSCs, which blocks the "positive inflammatory loop" in an IL-1Ra-dependent manner, and induces macrophage polarization toward the M2 phenotype [34][35][36][37] . This is the theoretical basis of stem cell-based therapy for arthritis. In this study, we explored the effects of IL-1β on lncRNA and mRNA expression profiles of human SFMSCs, which may help us understand in more detail the regulation mechanism underlying this biological process. www.nature.com/scientificreports/ The amount of differentially expressed lncRNAs and the lncRNA profile signature are likely affected by the exposure time of SFMSCs to inflammatory conditions. Here, we chose 2 h and not days because this exposure protocol was usually performed in previous study to explore the molecular function of lncRNAs 38 . In addition, we used the same method in our previous study to explored the molecular mechanism of how IL-1β affects biological behaviors of SFMSCs 6 . To further explore the underlying molecular mechanism, we used the same exposure method in this study.
In our previous study 6 , we also performed a dose-response experiment of IL-1β (0, 0.1, 1, and 10 ng/mL) and found inflammatory activation and SFMSCs concentration were positively correlated. In addition, the chondrogenic differentiation inhibition of SFMSCs by IL-1β was the strongest at 10 ng/mL. Therefore, we chose this concentration of IL-1β as it may be more helpful for us to explore the underlying mechanism of chondrogenic differentiation inhibition of SFMSCs by IL-1β.
To clarify the core subnetwork, signaling pathways, and core proteins involved in this biological process that SFMSCs in response to IL-1β, we selected differentially expressed mRNAs and performed GO and KEGG pathway analysis. We found that the differentially expressed mRNAs were mainly upregulated according to the results of volcano plots and heatmaps. Moreover, differentially expressed mRNAs showing significant enrichment in GO terms or KEGG pathways were upregulated, and few downregulated genes showed significant enrichment in GO terms or KEGG pathways, indicating that SFMSCs are activated by IL-1β stimulation. In addition, we analyzed the most significantly enriched GO terms in the biological process category and found that GO terms were mainly related to inflammatory activation, consistent with our in vitro analyses.
We also evaluated the distribution of the differentially expressed mRNAs with the top 30 K-scores and found that these mRNAs were enriched in chromosome 4, suggesting that IL-1β may have the most substantial effects on genes found on chromosome 4. Interestingly, based on analysis of the GO terms belonging to the cellular component category, we found that the differentially expressed mRNAs were only significantly enriched in the GO term IκB/NF-κB complex. In addition, from KEGG pathway analysis, we found that NF-κB pathway was one of the most significantly enriched pathways. Protein-protein interaction network analysis using the STRING database showed that some protein components of the NF-κB pathway were involved in the subnetwork with the highest K-score. Altogether, these results suggest that NF-κB pathway activation may be a core pathway in this biological process, consistent with our previous studies 6, 39 .
LncRNAs play important regulatory roles in many biological processes. In this study, our microarray data identified many differentially expressed lncRNAs, indicating that lncRNAs are involved in IL-1β stimulation of SFMSCs. However, the functions of most lncRNAs remain unclear. Analyzing the relationships between lncRNAs and paired neighboring mRNAs is a common method for annotating lncRNAs. In this study, antisense lncRNAs, lincRNAs, or enhancer lncRNAs and paired mRNAs that showed differential expression were selected. We then analyzed one pair containing a downstream molecule from the NF-κB pathway (LOC541472, IL-6) and found that the lncRNA LOC541472 was closely related to IL-6; both genes were significantly upregulated after IL-1β stimulation, and knockdown of LOC541472 also downregulated IL-6 expression and blocked NF-κB pathway activation, indicating that LOC541472 may be closely involved in regulating IL-6 expression and the NF-κB pathway. RNAscope analysis indicated that LOC541472 may be related to IL-6 expression and its translocation behavior from the nucleus into the cytoplasm. Further analyses of the mechanisms mediating the interactions between these molecules are required.
The functions of lncRNAs are commonly annotated by analyzing the functions of mRNAs that are coexpressed with the lncRNA. To provide some insights into the functions of lncRNAs in SFMSCs after IL-1β stimulation, we performed weighted correlation network analysis for the differentially expressed lncRNAs and mRNAs in SFMSCs after IL-1β stimulation. We visualized the highest K-score subnetwork, which contained three lncRNAs and three mRNA. Interestingly, this subnetwork contained one of the NF-κB pathway downstream molecules (IL-8). We then selected the lncRNA with the highest K-score (RP3-467K16.4) in this subnetwork and verified its function by analysis of GO terms and KEGG pathways associated with co-expressed mRNAs. We found significant enrichment of GO terms related to inflammatory activation and KEGG pathways related to the NF-κB pathway. Further evaluations showed that RP3-467K16.4 and IL-8 were both significantly upregulated after IL-1β stimulation, and knockdown of RP3-467K16.4 downregulated IL-8 expression and blocked NF-κB pathway activation. Additional studies are required to further assess the relationships among RP3-467K16.4, IL-8, and the NF-κB signaling pathway.
Both immunoregulatory function and cartilage differentiation potential of SFMSCs are important repair mechanisms for temporomandibular joint damage and joint homeostasis. In this study, we only investigated the lncRNA and mRNA profiles in untreated and IL-1β-inflamed SFMSCs. Our findings may be helpful for investigating the effect of lncRNAs and mRNAs on SFMSCs in this process. However, our study is limited in terms of investigating the immunoregulatory function of SFMSCs. As is known, MSCs can affect inflammation-related players by secreting cytokines, such as IL6 and IL8, for immunoregulation. In addition, secretion of exosomes also represents an important pathway for MSCs to participate in immune regulation. It is known that exosomes can be internalized by macrophages-the major inflammation-related player-and that the RNAs contained in exosomes can regulate macrophage polarization [40][41][42] . Therefore, further research should investigate the exosomal RNAs profile of IL-1β untreated/inflamed SFMSCs.

Conclusion
There were significant differences in lncRNA and mRNA expression profiles after treatment of SFMSCs with IL-1β, and the NF-κB pathway may be one of the most crucial pathways in this biological process. We found a co-expression relationship between RP3-467K16. 4 and IL8 and between LOC541472 and IL6, which are related