Aberrant lncRNA–mRNA expression profile and function networks during the adipogenesis of mesenchymal stem cells from patients with ankylosing spondylitis

Objective: We have already demonstrated that mesenchymal stem cells from patients with ankylosing spondylitis (ASMSCs) exhibited greater adipogenic differentiation potential than those from healthy donors (HDMSCs). Here, we further investigated the expression profile of long noncoding RNA (lncRNA) and mRNA, aiming to explore the underlying mechanism of abnormal adipogenic differentiation in ASMSCs.

Methods: HDMSCs and ASMSCs were separately isolated and induced with adipogenic differentiation medium for 10 days. Thereafter, lncRNAs and mRNAs that were differentially expressed (DE) between HDMSCs and ASMSCs were identified via high-throughput sequencing and confirmed by quantitative real-time PCR (qRT–PCR) assays. Then, the DE genes were annotated and enriched by GO analysis. In addition, protein interaction network was constructed to evaluate the interactions between DE mRNAs and to find hub nodes and study cliques. Besides, co-expression network analysis was carried out to assess the co-expressions between DE mRNA and DE lncRNAs, and competing endogenous RNA (ceRNA) network analysis were conducted to predict the relationships among lncRNAs, mRNAs and miRNAs. The signaling pathways based on the DE genes and the predicted DE genes were enriched by KEGG analysis.

Results: A total of 263 DE lncRNAs and 1376 DE mRNAs were found during adipogenesis in ASMSCs. qRT–PCR indicated that the expression of the top 20 mRNAs and the top 10 lncRNAs was consistent with the high-throughput sequencing data. Several lncRNAs (NR_125386.1, NR_046473.1 and NR_038937.1) and their target genes (SPN and OR1AIP2), together with the significantly co-expressed pairs of DE lncRNAs and DE mRNAs (SLC38A5-ENST00000429588.1, TMEM61-ENST00000400755.3 and C5orf46-ENST00000512300.1), were closely related to the enhanced adipogenesis of ASMSCs by modulating the PPAR signaling pathway.

Conclusion: Our study analyzed the expression profiles of DE lncRNAs and DE mRNAs during adipogenesis in ASMSCs and HDMSCs. Several DE lncRNAs, DE mRNAs and signaling pathways that probably participate in the aberrant adipogenesis of ASMSCs were selected for future study. These results will likely provide potential targets for our intervention on fat metaplasia and subsequent new bone formation in patients with AS in the future.

Introduction

Ankylosing spondylitis (AS) is a common inflammatory arthritis characterized by chronic inflammation and pathological new bone formation (Navarro-Compan et al., 2021). With the use of anti-tumor necrosis factor-α agents, nonsteroidal anti-inflammatory drugs and disease-modifying anti-rheumatic drugs, low back pain and peripheral joint pain caused by chronic inflammation can be significantly relieved in most patients (Ward et al., 2019). However, structural damage and disability, resulting from pathological new bone formation, seem inevitable for them (Ritchlin and Adamopoulos, 2021). Thus, further study the pathogenesis of new bone formation and provide new treatment options are urgently needed.

Fat metaplasia, with enhanced signals on T1-weighted sequences and reduced signals on short tau inversion recovery sequences, is an important magnetic resonance imaging (MRI) finding in patients with AS (Maksymowych et al., 2014). It can be frequently observed in the sacroiliac joint or the vertebral corner, where inflammation and new bone formation mainly occur (Machado et al., 2016). Recently, evidence has shown that this MRI signal corresponds to a high level of adipocyte accumulation (Baraliakos et al., 2019), which is a vital intermediary step and a strong contributor to pathological new bone formation in AS (Chiowchanwisawakit et al., 2011). Therefore, elucidating the cause of fat metaplasia is important.

Mesenchymal stem cells (MSCs) are a heterogeneous population of cells with fibroblast-like morphology and immunomodulation potentials, as well as self-renewal and multipotent differentiation capacities (Uccelli et al., 2008). As an origin of both fat and bone, MSCs obviously contribute to the abnormal adipocyte accumulation and thus the pathologic structural changes in AS. Our previous study demonstrated that MSCs from patients with AS (ASMSCs) had a significantly enhanced adipogenic capacity compared with those from healthy donors (HDMSCs), potentially leading to adipocyte accumulation and fat metaplasia (Liu et al., 2019). However, the specific mechanism of abnormal ASMSC adipogenesis remains undefined.

Long noncoding RNAs (lncRNAs) are RNA transcripts containing more than 200 nucleotides but without protein-coding potential (Ransohoff et al., 2018). Through interacting with DNA, RNAs and/or proteins, they play crucial roles in various cellular processes, including cell differentiation, the cell cycle and metabolism (Bridges et al., 2021). Evidence is accumulating that lncRNAs not only participate in the adipogenic differentiation of MSCs, but also contribute to the pathogenesis of AS (Sun et al., 2022). Based on these findings, it is reasonable for us to presume that the enhanced adipogenesis of ASMSCs may be related to abnormal lncRNA expression.

In this study, we compared the lncRNA and mRNA expression profiles of both ASMSCs and HDMSCs and analyzed the possible mechanism underlying the expression difference, aiming to reveal the role lncRNAs may play in the abnormal adipogenesis of ASMSCs and provide insights into the pathogenesis of new bone formation.

Materials and methods

Ethics and enrollment

Our study was approved by the Ethics Committee of Zhujiang Hospital, Southern Medical University, Guangzhou, China. After obtaining written informed consent, a punch in the posterior upper iliac crest was generated after local anesthesia with 2% lidocaine, and bone marrow samples were extracted from 12 healthy donors and 10 patients with AS (diagnosed according to the New York Revised Criteria (van der Linden et al., 1984)). Detailed characteristics of these volunteers are presented in Supplementary Table S1.

Cell isolation and culture

HDMSCs and ASMSCs were separately purified from the bone marrow samples of healthy donors and AS patients using density gradient centrifugation as described (Colter et al., 2000). To be specific, 20 ml bone marrow sample was diluted 1:1 with Hanks’ balanced salt solution (HBSS, GIBCO, United States) and layered over about 10 ml of Ficoll (Ficoll-Paque, Pharmacia, Sweden). After centrifugation (2,500 × g, 30 min), the mononuclear cell layer was extracted from the interface, suspended in HBSS, then centrifuged (1,500 × g, 15 min) and resuspended in Dulbecco’s modified Eagle’s medium (DMEM, Gibco, United States) supplemented with 10% heat-inactivated fetal bovine serum (FBS, Gibco, United States). Afterwards, they were seeded in 25-cm² flasks and cultured in an incubator with 5% CO₂ at 37°C. To remove the suspended cells, the medium was replaced 48 h later and every 3 days thereafter. At 80–90% confluence, the adherent cells were harvested using 0.25% Trypsin containing 0.53 mM EDTA (Gibco, United States) and re-plated in 75-cm² flasks as passage 1 of MSCs. MSCs were expanded and used for adipogenic differentiation at passages 4. Flow cytometry was performed to identify the surface markers of MSCs by detecting the expression of CD29-PE, CD14-APC, CD44-FITC, CD105-FITC, CD45-APC and HLA DR-PerCP (all the antibodies were purchased from BD, United States).

Adipogenic differentiation induction

HDMSCs and ASMSCs were seeded in 12-well plates at a density of 5 × 10⁴ cells/well in growth medium containing high-glucose DMEM supplemented with 10% FBS. When the cells reached 80–90% confluence, the medium was replaced with adipogenic differentiation medium consisting of high-glucose DMEM supplemented with 10% FBS, 0.2 mM indomethacin (Sigma–Aldrich, Germany), 10 μg/ml insulin (Sigma–Aldrich, Germany), 0.5 mM 3-isobutyl-1-methylxanthine (Sigma–Aldrich, Germany) and 1 μM dexamethasone (Sigma–Aldrich, Germany). The adipogenic differentiation medium was replaced every 3 days. HDMSCs and ASMSCs on the 10th day of induction were harvested for the following experiments.

Oil red O staining and Triglyceride quantification

As described in our previous study (Liu et al., 2019; Cen et al., 2020), cells were firstly washed with phosphate-buffered saline and then fixed with 4% formaldehyde for 20 min. Thereafter, formaldehyde was removed, and the cells were washed with 60% isopropanol and then stained with 0.2% ORO (Sigma-Aldrich, Germany) for 30 min. Afterwards, cells were washed with phosphate-buffered saline and observed under a microscope (Olympus, Japan). Stained oil droplets were dissolved in 100% isopropanol and quantified by detecting the optical absorbance (500 nm) with a spectrophotometer (Thermo Fisher, Germany). Triglyceride content was quantified using a triglyceride quantification kit (ab178780, Abcam, United Kingdom), according to manufacturer’s instructions.

Western blot analysis

Proteins were extracted from HDMSCs and ASMSCs on the 10th day of induction, and then quantified as previously described (Liu et al., 2019). Equal amounts of protein extracts were denatured, separated and then transferred onto polyvinylidene fluoride membranes (Millipore, United States). Afterwards, membranes were blocked with 5% skim milk and incubated overnight at 4°C with primary antibodies against GAPDH, peroxisome proliferator-activated receptor gamma (PPAR-γ), fatty acid binding protein 4 (FABP4), adiponectin (all from Abcam, United Kingdom). Horseradish peroxidase -conjugated immunoglobulin IgG (Santa Cruz, United States) was used as secondary antibody and incubated with the membranes at room temperature. The immunoreactive bands were then visualized using the Immobilon Western Chemiluminescent HRP Substrate (Millipore, United States).

Library construction and sequencing

On the 10th day of induction, total RNAs were respectively extracted from HDMSCs and ASMSCs using the TRIzol (Invitrogen, United States) according to the manufacturer’s instructions, and ribosomal RNA was removed using the Ribo-Zero™ kit (Epicentre, Madison, WI, United States). Fragmented RNA (the average length was approximately 200 bp) were subjected to first strand and second strand cDNA synthesis following by adaptor ligation and enrichment with a low-cycle according to instructions of NEBNext^® Ultra™ RNA Library Prep Kit for Illumina (NEB, United States). The purified library products were evaluated using the Agilent 2200 TapeStation and Qubit^®2.0(Life Technologies, United States). The libraries were paired-end sequenced (PE150, Sequencing reads were 150 bp) at Guangzhou RiboBio Co., Ltd. (Guangzhou, China) using Illumina-HiSeq 3000 platform.

Quality control and expression analysis

In order to remove the low-quality data, the Raw fastq sequences were treated with Trimmomatic tools (v 0.36) using the following options: TRAILING:20, SLIDINGWINDOW:4:15 MINLEN:52, to remove trailing sequences below a phred quality score of 20 and to achieve uniform sequence lengths for downstream clustering processes. Thereafter, the relative high-quality data were normalized to the expected number of reads per kilobase of transcript per million mapped reads (RPKM). Differentially expressed (DE) genes were detected by the DESeq2 method based on negative binomial generalized linear models. Genes with significant fold changes (log2-fold change >1; Q-value < 0.05) between HDMSCs and ASMSCs were accepted for further study.

Quantitative real-time polymerase chain reaction

To confirm the reliability of RNA sequencing, qRT–PCR was performed as previously described (Liu et al., 2017). Simply put, total RNAs were extracted from HDMSCs and ASMSCs on the 10th day of induction, and then transcribed into complementary DNA by PrimeScript RT reagent kit (Takara, Japan). qRT–PCR was conducted with a Light Cycler^® 480 PCR System (Roche, Switzerland). Data were normalized to GAPDH, and the relative expression level of each gene was analyzed using the 2^−ΔCt method.

GO and KEGG pathway analysis

The annotation and enrichment information of DE genes were analyzed by Gene Ontology (GO) term enrichment analysis using KOBAS3.0 software (http://www.geneontology.org). Label classification for gene function and gene product attributes, including cellular component (CC), molecular function (MF) and biological process (BP), were detailed. The biological pathways of the DE genes were categorized by KEGG pathway analysis using KOBAS 3.0 software (http://www.genome.jp/kegg). After being calculated using Fisher’s exact test, signal transduction and disease pathway enrichment of DE genes with p values <0.05 were mapped using KEGG pathway annotation.

Interaction analysis and co-expression network construction

To obtain information on all target-related genes or proteins, PPI network was constructed by using the STRING database and visualized in Cytoscape software as described (confidence score was set as score >0.4) (Doncheva et al., 2019). GENEMANIA (http://genemania.org/search/) was used to evaluate the interactions between them (Warde-Farley et al., 2010). Additionally, the MCODE app was used to find hub nodes and study cliques in the PPI network (degree cutoff = 2, max. Depth = 100, k-core = 2, and node score cutoff = 0.2) (Bader and Hogue, 2003). After obtaining the overall interaction relationship, the Pearson correlation coefficients and p values between multiple genes were calculated. lncRNAs and mRNAs with COR ≤0.85 and p value ≥ 0.05 were excluded, and the mRNA–lncRNA co-expression network was also built by Cytoscape software (https://cytoscape.org).

Competing endogenous RNA network analysis

The above-selected mRNAs and lncRNAs for co-expression network analysis were also used to predict miRNA in the miRbase. Then, miRNAs obtained from miRbase were further screened using the miRanda and TargetScan programs. Afterward, RNA22 was used to predict those lncRNAs and mRNAs with miRNA recognition elements (MREs) for the targeted miRNAs. Finally, we constructed the competitive ceRNA network using Cytoscape software.

Statistical analysis

Quantitative data are expressed as the means ± standard deviations (SD). Spearman correlation was used to assess the relationship between lncRNAs and their target genes. Student’s t test was performed to analyze the statistical significance between two groups. The abovementioned statistical analysis was performed using SPSS 23.0 software (Chicago), and p values < 0.05 indicated statistical significant.

Results

ASMSCs had greater adipogenic capacity than HDMSCs

Both HDMSCs and ASMSCs used in our study were plastic-adherent cells with spindle-shaped, and they were negative for CD14, CD45 and HLA-DR and positive for CD29, CD44 and CD105 (Supplementary Figure S1), which conformed to the criteria stated by the International Society for Cellular Therapy (Dominici et al., 2006). HDMSCs and ASMSCs were cultured in adipogenic medium for 0–10 days. Cells were stained with ORO on Days 0 and 10, and the stained oil droplets were dissolved in 100% isopropanol. ASMSCs exhibited more intense staining on Day 10 than HDMSCs (Figure 1A). Consistent results were also observed in the quantification of ORO staining (Figure 1B) and triglyceride accumulation (Figure 1C). Moreover, protein expressions of the adipogenic differentiation markers (PPAR-γ, FABP4, adiponectin) were higher in ASMSCs than those of HDMSCs (Figure 1D). As with our previous report (Liu et al., 2019), these data supported the notion that ASMSCs had greater adipogenic capacity than HDMSCs.

FIGURE 1

FIGURE 1. ASMSCs had greater adipogenic capacity than HDMSCs. The adipogenic differentiation potential of HDMSCs and ASMSCs were assessed via ORO staining and quantification and further confirmed by quantification of triglyceride accumulation. (A) ORO staining indicated that ASMSCs displayed more intense staining on Day 10 (100X, scale bar indicate 100 um). Consistent results were also observed in quantifying of ORO staining (B) and triglyceride accumulation (C). Moreover, protein expressions of PPAR-γ, FABP4 and adiponectin were higher in ASMSCs than those of HDMSCs (D). *p < 0.05.

Expression profile of DE mRNAs and DE lncRNAs during adipogenesis

As depicted in the clustergram (Figure 2A) and volcano plot (Figure 2B), a total of 1376 DE mRNAs were found between HDMSCs and ASMSCs. Among them, 630 mRNAs were upregulated, while 746 mRNAs were downregulated in ASMSCs. As previously described (Li M. et al., 2019), the top 20 DE mRNAs with largest fold change were chosen for further analysis and listed in Table 1. Adipogenesis-related mRNAs, including GJB2, ADRA1A, SYT1, LPO, BICDL2, TBATA, APELA, CCKAR, RAPGEF3, NTSR1 and CDH10, may be related to abnormal adipogenesis of ASMSCs. As shown in the clustergram (Figure 2C) and volcano plot (Figure 2D), 137 upregulated and 126 downregulated lncRNAs were also found in ASMSCs. As previously described (Wang et al., 2021), the top 10 DE lncRNAs with largest fold change were chosen for further analysis and shown in Table 2. These data revealed a significant genetic difference between HDMSCs and ASMSCs during adipogenesis.

FIGURE 2

FIGURE 2. Expression profile of DE mRNAs and DE lncRNAs during adipogenesis. (A) Heatmaps of DE mRNAs between HDMSCs and ASMSCs. (B) Volcano plots of DE mRNAs between HDMSCs and ASMSCs. A total of 630 upregulated mRNAs and 746 downregulated mRNAs were found. (C) Heatmaps of DE lncRNAs between HDMSCs and ASMSCs. (D) Volcano plots of DE lncRNAs between HDMSCs and ASMSCs. A total of 137 upregulated mRNAs and 126 downregulated lncRNAs were found. In the heatmap, HD indicates MSCs from healthy donors, and AS indicates MSCs from patients with ankylosing spondylitis.

TABLE 1

TABLE 1. The characteristics of the top 20 mRNA with largest fold change.

TABLE 2

TABLE 2. The characteristics of the top 10 LncRNA with largest fold change.

Confirmation of DE mRNAs and lncRNAs by qRT–PCR

To detect the reliability of the RNA sequencing, the expression levels of the top 20 DE mRNAs and the top 10 DE lncRNAs were detected by qRT–PCR. Compared with HDMSCs, ADRA1A, TBATA, CCDC144NL, NXNL1, CFAP47 and CDH10 were upregulated in ASMSCs (Figure 3A), while GJB2, KLHL4, SYT1, etc., were downregulated in ASMSCs (Figure 3B). With regard to lncRNAs, the expression levels of ENST100000546836.1, NR_023925.1, ENST00000432265.1 and NR_051996.1 were increased in ASMSCs (Figure 3C), while the expression levels of ENST00000593060.1, ENST00000592405.1, ENST00000429588.1, ENST00000400755.3, ENST00000512300.1 and NR_003948.2 were decreased in ASMSCs (Figure 3D). These results were essentially comparable to the data shown in Table 1 and Table 2, which strongly confirmed the reliability and accuracy of the RNA sequence data.

FIGURE 3

FIGURE 3. Confirmation of DE mRNAs and DE lncRNAs by qRT–PCR. The expression levels of the top 20 DE mRNAs and top 10 DE lncRNAs were confirmed via qRT–PCR. As shown in (A) and (B), 6 mRNAs were upregulated and 14 mRNAs were downregulated in ASMSCs. As depicted in (C) and (D), 4 lncRNAs were upregulated and 6 lncRNAs were downregulated in ASMSCs. Data are presented as the means ± SDs. *p < 0.05.

GO and KEGG pathway analysis

DE mRNAs between HDMSCs and ASMSCs were annotated and enriched by GO analysis. The top 10 GO terms of the three domains, including biological processes (BP), cellular components (CC) and molecular function (MF), are presented in Figure 4A and Table 3. Specifically, in the BP domain, the significantly enriched GO terms included multicellular organism development, anatomical structure morphogenesis, organic acid metabolic process, etc. In the CC domain, the significantly enriched GO terms were collagen-containing extracellular matrix, plasma membrane, lipid droplet, etc. In the MF domain, the significantly enriched GO terms were protein homodimerization activity, rho guanyl-nucleotide exchange factor activity, ras guanyl-nucleotide exchange factor activity, etc.

FIGURE 4

FIGURE 4. GO and KEGG pathway analysis. (A) Top 10 GO terms with the largest significant difference in the three domains. (B) Top 10 KEGG signaling pathways with the largest significant difference.

TABLE 3

TABLE 3. Go analysis of the top 10 mRNA expression with largest significant difference in three domains.

KEGG pathway analysis was performed to enrich the critical signaling pathways involved in the abnormal adipogenesis of ASMSCs. A total of 381 signaling pathways were enriched, and the top 10 pathways with the largest significant difference are shown in Figure 4B and Table 4. They are fatty acid metabolism, amino acid biosynthesis, carbon metabolism, peroxisome, alanine aspartate and glutamate metabolism, glycolysis/gluconeogenesis, regulation of lipolysis in adipocytes, PPAR signaling pathway, AMPK signaling pathway, and HIF-1 signaling pathway. Among them, the PPAR signaling pathway and related mRNAs might play important roles in the increased adipogenesis of ASMSCs.

TABLE 4

TABLE 4. The top 10 pathways with largest significant difference in KEGG analysis.

Interaction and co-expression network analysis

Illustrating the molecular interaction helps to clarify the pathophysiology and reveal the underlying mechanisms of disease. Proteins encoded by DE mRNAs were analyzed using PPI network. Our data suggested that GRIA1, FLT1, APOE, ENO2, GNG4, NTRK3 and SCN1A interacted with much more DE genes than others, implying that they might be important targets contributing to the abnormal adipogenesis of AS (Figure 5A and Supplementary Table S2). Meanwhile, the specific interactions between DE genes were shown in Figure 5B and Supplementary Table S3. A total of 2349 interactions between DE genes were identified. Among them, 1275 were co-expression, 732 were genetic interaction, 111 were physical interaction, 106 were shared protein domains, 89 were co-localization, 34 were pathway-related and 2 were predicted. Besides, module analysis was also performed to find and study cliques in the PPI network. As shown in Figure 5C and Supplementary Table S4, the top 3 clusters and their hub nodes were: GJA3, GJB2 and GJA5 cluster, AREG, TGFA and FLT1 cluster, ENO2, ADH1A, APOC1, ADH1B, APOE and ABCG1 cluster.

FIGURE 5

FIGURE 5. Interaction and co-expression network analysis. (A) Protein–protein interaction (PPI) network based on DE mRNAs. (B) The specific interactions between DE genes in the PPI network. (C) The significant hub nodes and molecular complex clustered by the MCODE in the PPI network. (D) Co-expression network of DE lncRNAs and DE mRNAs. (E) ceRNA network of DE lncRNAs, mRNAs and their predicted miRNAs.

Thereafter, interactions between DE lncRNAs and DE mRNAs were analyzed, and the top 30 pairs are shown in Figure 5D after filtering the repetitions. Within the network, SLC7A5, SLC38A, FGF11, TMEM61, NPAS2, HIBCH, c21orf33, PNPLA3, FCRLA and their co-expressed lncRNAs were significantly enriched. The top 10 co-expression pairs are listed in Table 5. ceRNA network analysis is necessary to reveal the interactions among lncRNAs, mRNAs and miRNAs. As shown in Figure 5E, 32 predicted miRNAs were significantly enriched and combined with 3 lncRNAs and 182 mRNAs. Among them, miR-6778-5p, miR-6127, miR-6089, miR-6813-5p and miR-149-3p interacted with much more mRNAs than others, indicating that they might play critical roles in the mechanism of enhanced adipogenesis in ASMSCs.

TABLE 5

TABLE 5. The top 10 co-expression pairs.

Target prediction of lncRNAs

lncRNAs regulate gene expression by interacting with proteins, RNAs and DNAs. The possible target genes of the DE lncRNAs were analyzed in this study, and the predictive genes with binding scores greater than 0.9 are shown in Figure 6A. NR_125386.1, NR_046473.1 and NR_038937.1 were the top 3 lncRNAs that appeared to have various target genes. SPN and OR1AIP2 were target genes of all the top 3 lncRNAs. HAUS2, TRPV1, PTCHD4, SKA1, NQO1, H6PD, ORAI2, RRP15, COX6B2, SLC35F6, FBXL18, PIGW, LPCAT2, FUT2, PNMA2, ZFP42, CBX5, UTP11 and TRAF3IP2 were target genes of two lncRNAs. In addition, the Venn diagram analysis showed that 25 genes existed both in 553 DE mRNAs and in 548 DE lncRNA target genes (Figure 6B). Moreover, the target genes of DE lncRNAs were clustered by KEGG pathway analysis (Figure 6C). The PPAR signaling pathway, adipocytokine signaling pathway, Wnt signaling pathway and the other 12 signaling pathways related to DE lncRNAs might play critical roles in the abnormal adipogenesis of ASMSCs.

FIGURE 6

FIGURE 6. Target prediction of lncRNAs. (A) DE lncRNAs and their target genes. Blue indicates downregulated, while orange indicates upregulated. (B) Intersection between DE mRNAs and target genes of DE lncRNAs. (C) KEGG pathway analysis based on the DE lncRNA-targeted genes.

Discussion

The increased adipogenesis of ASMSCs may result in fat metaplasia and contribute to new bone formation. Further investigation of the underlying mechanism of increased adipogenesis of ASMSCs may help to elucidate the pathogenesis of AS. In this study, we performed high-throughput sequencing to analyze lncRNA and mRNA expression in ASMSCs and HDMSCs during adipogenesis. Our data demonstrated that a total of 263 lncRNAs and 1376 mRNAs were abnormally expressed, and the PPAR signaling pathway and its related mRNAs might be involved in the increased adipogenesis of ASMSCs. Further bioinformatics analysis indicated that several DE lncRNAs (ENST00000429588.1, ENST00000400755.3 & ENST00000512300.1) may be upstream targets contributing to abnormal ASMSC adipogenesis by acting as ceRNAs or by interacting with mRNAs.

MSCs are recognized as a promising and revolutionary treatment due to their multiple differentiation potential and immunomodulatory properties. An increasing number of studies have reported that MSCs can effectively alleviate inflammation and other symptoms of rheumatism (Shi et al., 2018). However, the abnormal behavior of MSCs may also result in autoimmune diseases, including AS (Krajewska-Wlodarczyk et al., 2017). Ye et al. reported that oxidative stress-mediated mitochondrial dysfunction facilitates senescence of ASMSCs and contributes to the pathogenesis of AS (Ye et al., 2020). Xie et al. demonstrated that TNF-α-mediated m⁶A modification of ELMO1 triggers directional migration of ASMSCs, which might lead to chronic inflammation of AS (Xie et al., 2021). Similarly, our previous data has shown that ASMSCs have enhanced adipogenic capacity (Liu et al., 2019). Given that an imbalance between adipogenesis and osteogenesis may lead to abnormal bone metabolism (Cen et al., 2020), our previous study may partly explain the phenomenon of fat metaplasia and new bone formation in patients with AS. Therefore, further study of the underlying mechanism of abnormal adipogenesis in ASMSCs is of significance.

lncRNAs have aroused great interest in medical and scientific circles due to their multifaceted and versatile regulatory roles in various biological processes. Numerous lncRNAs are reported to modulate the adipogenic differentiation process (Squillaro et al., 2020). For example, SRA, lnc-ORA and lncRNA-Adi can regulate the early stage of adipogenesis, while ADNCR, lncRNA Plnc1, lncRNA PVT1 and NEAT1 were suggested to target C/EBPα and/or PPARγ in the late stage of adipogenesis (Rey et al., 2021). In addition, lncRNAs have also been reported to be involved in the pathogenesis of AS and serve as biomarkers and/or potential therapeutic targets (Chen et al., 2019). For example, lncRNA-AK001085 and TUG1 can serve as diagnostic biomarkers for AS. LINC00311, NKILA and Lnc-ITSN1-2 can be used to monitor the activity and assess the prognosis of AS. H19 and MEG3 are potential therapeutic targets for AS (Sun et al., 2022). In our study, 137 upregulated and 126 downregulated lncRNAs were found in ASMSCs during adipogenesis as compared to HDMSCs. Among them, three of the top 10 DE lncRNAs (ENST00000429588.1, ENST00000400755.3 & ENST00000512300.1) appeared to be significantly co-expressed with DE mRNAs, implying that they may be important targets of abnormal adipogenesis. However, their specific role still needs to be further confirmed.

Considering that lncRNAs play regulatory roles in gene expression by interacting with various molecular species, bioinformatic analysis was performed to illuminate their interactions and provide directions for our future study. The PPI network analysis indicated that, in addition to our previously confirmed genes (PPARG, CEBPA and ADIPOQ) (Liu et al., 2019), other DE genes (GRIA1, FLT1, APOE, ENO2, GNG4, NTRK3 and SCN1A) warrant further study as well. Among them, it is worth noting that APOE was also a hub node gene in the modules analysis. Given the fact that APOE is a key player in adipogenesis (Chiba et al., 2003), it is reasonable to presume that APOE is involved in the abnormal adipogenesis of ASMSCs. Additionally, the co-expression analysis provided three promising co-expression candidates (SLC38A5-ENST00000429588.1, TMEM61-ENST00000400755.3 and C5orf46-ENST00000512300.1) involved in abnormal adipogenesis. Moreover, the ceRNA network analysis suggested that miR-6778-5p, miR-6127, miR-6089, miR-6813-5p and miR-149-3p were significantly enriched. Notably, miR-149-3p can regulate both the adipogenic and osteogenic differentiation of MSCs(Li Y. et al., 2019), while miR-6127 and miR-6089 are related to immunoregulation and inflammation (Tubita et al., 2019; Yang et al., 2020). Finally, target prediction analysis revealed that the lncRNAs (NR_125386.1, NR_046473.1 and NR_038937.1) and the target genes they shared (SPN and OR1AIP2) may play critical roles in the abnormal adipogenesis of ASMSCs through the PPAR signaling pathway.

The PPAR signaling pathway is one of the most important transcriptional modulators involved in both the early and late stages of adipogenic differentiation (Lee et al., 2019). In this study, the PPAR signaling pathway and its related mRNAs also drew our attention. KEGG analysis based on both DE mRNA- and DE lncRNA-targeted genes demonstrated that the PPAR signaling pathway was aberrantly activated and significantly enriched during adipogenesis in ASMSCs. Another study also found that serum from patients with AS could increase the expression of PPAR, the key molecule of the PPAR signaling pathway (Hu et al., 2015). Based on the recognized role of PPAR in adipogenic differentiation (Lefterova et al., 2014), we believe that the PPAR signaling pathway is an important cause and key target of fat metaplasia in AS. In addition, fatty acid metabolism (Xu et al., 2015), the HIF-1 signaling pathway (Ding et al., 2018) and alanine aspartate and glutamate metabolism (Zhou et al., 2020) were also reported to be related to AS, but further studies are still needed to explore whether they mediate abnormal adipogenesis.

In this study, we again confirmed the enhanced adipogenesis of ASMSCs. During adipogenesis, hundreds of DE lncRNAs and DE mRNAs were found in ASMSCs, and the potential regulatory mechanisms were explored by bioinformatics analysis. Several lncRNAs (NR_125386.1, NR_046473.1 and NR_038937.1) and their target genes (SPN and OR1AIP2), together with several co-expression pairs of DE lncRNAs and DE mRNAs (SLC38A5-ENST00000429588.1, TMEM61-ENST00000400755.3 and C5orf46-ENST00000512300.1), may play a central role in regulating the adipogenic differentiation of ASMSCs by activating the PPAR signaling pathway. These results could help to elucidate the pathogenesis of fat metaplasia and new bone formation in patients with AS and provide potential therapeutic targets for them. However, there are still some limitations in this study. For example, we did not clarify the specific function and relative mechanism of DE lncRNAs during adipogenesis. These limitations should be explored in depth in future studies.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: Gene Expression Omnibus accession number GSE212613. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE212613

Ethics statement

The studies involving human participants were reviewed and approved by The Ethics Committee of Zhujiang Hospital, Southern Medical University. Written informed consent to participate in this study was provided by the participants’ legal guardian/next of kin.

Author contributions

ZL and SQ designed this study, analyzed the data and revised the manuscript, they are the corresponding authors. SC, MC, and YW carried out this study and wrote the manuscript, they contributed equally to this work and share first authorship. MSCs were cultured and induced by XL and ZC. HC, YF, and CW help to perform the PCR and analyzed part of the data.

Funding

This study was financially supported by the Natural Science Foundation of Guangdong Province (2018A0303130258, 2022A1515012021 & 2020A1515110930) and the Scientific Research Project of Southern Medical University in China (PY2018N056).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.991875/full#supplementary-material

Supplementary Figure S1 | HDMSCs and ASMSCs exhibited similar morphologies and phenotypes. (A) Both HDMSCs and ASMSCs were plastic-adherent cells with spindle-shaped. (B) HDMSCs and ASMSCs were both negative for CD14, CD45 and HLA-DR and positive for CD29, CD44 and CD105, indicating a typical MSC phenotype.

References

Bader, G. D., and Hogue, C. W. (2003). An automated method for finding molecular complexes in large protein interaction networks. Bmc Bioinforma. 4, 2. doi:10.1186/1471-2105-4-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Baraliakos, X., Boehm, H., Bahrami, R., Samir, A., Schett, G., Luber, M., et al. (2019). What constitutes the fat signal detected by mri in the spine of patients with ankylosing spondylitis? A prospective study based on biopsies obtained during planned spinal osteotomy to correct hyperkyphosis or spinal stenosis. Ann. Rheum. Dis. 78, 1220–1225. doi:10.1136/annrheumdis-2018-214983

PubMed Abstract | CrossRef Full Text | Google Scholar

Bridges, M. C., Daulagala, A. C., and Kourtidis, A. (2021). Lnccation: Lncrna localization and function. J. Cell Biol. 220, e202009045. doi:10.1083/jcb.202009045

PubMed Abstract | CrossRef Full Text | Google Scholar

Cen, S., Li, J., Cai, Z., Pan, Y., Sun, Z., Li, Z., et al. (2020). Traf4 acts as A fate checkpoint to regulate the adipogenic differentiation of mscs by activating Pkm2. Ebiomedicine 54, 102722. doi:10.1016/j.ebiom.2020.102722

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, C., Rong, T., Li, Z., and Shen, J. (2019). Noncoding rnas involved in the pathogenesis of ankylosing spondylitis. Biomed. Res. Int. 2019, 6920281. doi:10.1155/2019/6920281

PubMed Abstract | CrossRef Full Text | Google Scholar

Chiba, T., Nakazawa, T., Yui, K., Kaneko, E., and Shimokado, K. (2003). Vldl induces adipocyte differentiation in apoe-dependent manner. Arterioscler. Thromb. Vasc. Biol. 23, 1423–1429. doi:10.1161/01.ATV.0000085040.58340.36

PubMed Abstract | CrossRef Full Text | Google Scholar

Chiowchanwisawakit, P., Lambert, R. G., Conner-Spady, B., and Maksymowych, W. P. (2011). Focal fat lesions at vertebral corners on magnetic resonance imaging predict the development of new syndesmophytes in ankylosing spondylitis. Arthritis Rheum. 63, 2215–2225. doi:10.1002/art.30393

PubMed Abstract | CrossRef Full Text | Google Scholar

Colter, D. C., Class, R., Digirolamo, C. M., and Prockop, D. J. (2000). Rapid expansion of recycling stem cells in cultures of plastic-adherent cells from human bone marrow. Proc. Natl. Acad. Sci. U. S. A. 97, 3213–3218. doi:10.1073/pnas.97.7.3213

PubMed Abstract | CrossRef Full Text | Google Scholar

Ding, M., Guan, T. J., Wei, C. Y., and Chen, B. H. (2018). Identification of pathways significantly associated with spondyloarthropathy/ankylosing spondylitis using the subpathway method. Mol. Med. Rep. 18, 3825–3833. doi:10.3892/mmr.2018.9395

PubMed Abstract | CrossRef Full Text | Google Scholar

Dominici, M., Le Blanc, K., Mueller, I., Slaper-Cortenbach, I., Marini, F., Krause, D., et al. (2006). Minimal criteria for defining multipotent mesenchymal stromal cells. The international society for cellular Therapy position statement. Cytotherapy 8, 315–317. doi:10.1080/14653240600855905

PubMed Abstract | CrossRef Full Text | Google Scholar

Doncheva, N. T., Morris, J. H., Gorodkin, J., and Jensen, L. J. (2019). Cytoscape stringapp: Network analysis and visualization of proteomics data. J. Proteome Res. 18, 623–632. doi:10.1021/acs.jproteome.8b00702

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu, Z., Lin, D., Qi, J., Qiu, M., Lv, Q., Li, Q., et al. (2015). Serum from patients with ankylosing spondylitis can increase ppard, fra-1, Mmp7, opg and rankl expression in Mg63 cells. Clin. (Sao Paulo) 70, 738–742. doi:10.6061/clinics/2015(11)04

PubMed Abstract | CrossRef Full Text | Google Scholar

Krajewska-Wlodarczyk, M., Owczarczyk-Saczonek, A., Placek, W., Osowski, A., Engelgardt, P., and Wojtkiewicz, J. (2017). Role of stem cells in pathophysiology and Therapy of spondyloarthropathies-new therapeutic possibilities? Int. J. Mol. Sci. 19, E80. doi:10.3390/ijms19010080

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, J. E., Schmidt, H., Lai, B., and Ge, K. (2019). Transcriptional and epigenomic regulation of adipogenesis. Mol. Cell. Biol. 39, e00601. doi:10.1128/MCB.00601-18

PubMed Abstract | CrossRef Full Text | Google Scholar

Lefterova, M. I., Haakonsson, A. K., Lazar, M. A., and Mandrup, S. (2014). PPARγ and the global map of adipogenesis and beyond. Trends Endocrinol. Metab. 25, 293–302. doi:10.1016/j.tem.2014.04.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, M., Xie, Z., Cai, Z., Su, F., Zheng, G., Li, J., et al. (2019a). Lncrna-mrna expression profiles and functional networks of mesenchymal stromal cells involved in monocyte regulation. Stem Cell Res. Ther. 10, 207. doi:10.1186/s13287-019-1306-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y., Yang, F., Gao, M., Gong, R., Jin, M., Liu, T., et al. (2019b). Mir-149-3p regulates the switch between adipogenic and osteogenic differentiation of bmscs by targeting fto. Mol. Ther. Nucleic Acids 17, 590–600. doi:10.1016/j.omtn.2019.06.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Z., Gao, L., Wang, P., Xie, Z., Cen, S., Li, Y., et al. (2017). Tnf-alpha induced the enhanced apoptosis of mesenchymal stem cells in ankylosing spondylitis by overexpressing trail-R2. Stem Cells Int. 2017, 4521324. doi:10.1155/2017/4521324

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Z., Wang, P., Cen, S., Gao, L., Xie, Z., Wu, X., et al. (2019). Increased Bmpr1a expression enhances the adipogenic differentiation of mesenchymal stem cells in patients with ankylosing spondylitis. Stem Cells Int. 2019, 4143167. doi:10.1155/2019/4143167

PubMed Abstract | CrossRef Full Text | Google Scholar

Machado, P. M., Baraliakos, X., Van Der Heijde, D., Braun, J., and Landewe, R. (2016). Mri vertebral corner inflammation followed by fat deposition is the strongest contributor to the development of new bone at the same vertebral corner: A multilevel longitudinal analysis in patients with ankylosing spondylitis. Ann. Rheum. Dis. 75, 1486–1493. doi:10.1136/annrheumdis-2015-208011

PubMed Abstract | CrossRef Full Text | Google Scholar

Maksymowych, W. P., Wichuk, S., Chiowchanwisawakit, P., Lambert, R. G., and Pedersen, S. J. (2014). Fat metaplasia and backfill are key intermediaries in the development of sacroiliac joint ankylosis in patients with ankylosing spondylitis. Arthritis Rheumatol. 66, 2958–2967. doi:10.1002/art.38792

PubMed Abstract | CrossRef Full Text | Google Scholar

Navarro-Compan, V., Sepriano, A., El-Zorkany, B., and Van Der Heijde, D. (2021). Axial spondyloarthritis. Ann. Rheum. Dis. 80, 1511–1521. doi:10.1136/annrheumdis-2021-221035

PubMed Abstract | CrossRef Full Text | Google Scholar

Ransohoff, J. D., Wei, Y., and Khavari, P. A. (2018). The functions and unique features of long intergenic non-coding rna. Nat. Rev. Mol. Cell Biol. 19, 143–157. doi:10.1038/nrm.2017.104

PubMed Abstract | CrossRef Full Text | Google Scholar

Rey, F., Urrata, V., Gilardini, L., Bertoli, S., Calcaterra, V., Zuccotti, G. V., et al. (2021). Role of long non-coding rnas in adipogenesis: State of the art and implications in obesity and obesity-associated diseases. Obes. Rev. 22, E13203. doi:10.1111/obr.13203

PubMed Abstract | CrossRef Full Text | Google Scholar

Ritchlin, C., and Adamopoulos, I. E. (2021). Axial spondyloarthritis: New advances in diagnosis and management. Bmj 372, M4447. doi:10.1136/bmj.m4447

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi, Y., Wang, Y., Li, Q., Liu, K., Hou, J., Shao, C., et al. (2018). Immunoregulatory mechanisms of mesenchymal stem and stromal cells in inflammatory diseases. Nat. Rev. Nephrol. 14, 493–507. doi:10.1038/s41581-018-0023-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Squillaro, T., Peluso, G., Galderisi, U., and Di Bernardo, G. (2020). Long non-coding rnas in regulation of adipogenesis and adipose tissue function. Elife 9, e59053. doi:10.7554/eLife.59053

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, R., Wang, X., Sun, X., Zhao, B., Zhang, X., Gong, X., et al. (2022). Emerging roles of long non-coding rnas in ankylosing spondylitis. Front. Immunol. 13, 790924. doi:10.3389/fimmu.2022.790924

PubMed Abstract | CrossRef Full Text | Google Scholar

Tubita, V., Segui-Barber, J., Lozano, J. J., Banon-Maneus, E., Rovira, J., Cucchiari, D., et al. (2019). Effect of immunosuppression in mirnas from extracellular vesicles of colorectal cancer and their influence on the pre-metastatic niche. Sci. Rep. 9, 11177. doi:10.1038/s41598-019-47581-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Uccelli, A., Moretta, L., and Pistoia, V. (2008). Mesenchymal stem cells in Health and disease. Nat. Rev. Immunol. 8, 726–736. doi:10.1038/nri2395

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Der Linden, S., Valkenburg, H. A., and Cats, A. (1984). Evaluation of diagnostic criteria for ankylosing spondylitis. A proposal for modification of the New York criteria. Arthritis Rheum. 27, 361–368. doi:10.1002/art.1780270401

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, S., Wang, Z., Su, H., Chen, F., Ma, M., Yu, W., et al. (2021). Effects of long-term culture on the biological characteristics and rna profiles of human bone-marrow-derived mesenchymal stem cells. Mol. Ther. Nucleic Acids 26, 557–574. doi:10.1016/j.omtn.2021.08.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Ward, M. M., Deodhar, A., Gensler, L. S., Dubreuil, M., Yu, D., Khan, M. A., et al. (2019). 2019 update of the American college of rheumatology/spondylitis association of America/spondyloarthritis research and treatment network recommendations for the treatment of ankylosing spondylitis and nonradiographic axial spondyloarthritis. Arthritis Care Res. 71, 1285–1299. doi:10.1002/acr.24025

CrossRef Full Text | Google Scholar

Warde-Farley, D., Donaldson, S. L., Comes, O., Zuberi, K., Badrawi, R., Chao, P., et al. (2010). The genemania prediction server: Biological network integration for gene prioritization and predicting gene function. Nucleic Acids Res. 38, W214–W220. doi:10.1093/nar/gkq537

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, Z., Yu, W., Zheng, G., Li, J., Cen, S., Ye, G., et al. (2021). Tnf-alpha-mediated M(6)A modification of Elmo1 triggers directional migration of mesenchymal stem cell in ankylosing spondylitis. Nat. Commun. 12, 5373. doi:10.1038/s41467-021-25710-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, W. D., Yang, X. Y., Li, D. H., Zheng, K. D., Qiu, P. C., Zhang, W., et al. (2015). Up-regulation of fatty acid oxidation in the ligament as A contributing factor of ankylosing spondylitis: A comparative proteomic study. J. Proteomics 113, 57–72. doi:10.1016/j.jprot.2014.09.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, J., Cheng, M., Gu, B., Wang, J., Yan, S., and Xu, D. (2020). CircRNA_09505 aggravates inflammation and joint damage in collagen-induced arthritis mice via miR-6089/AKT1/NF-κB axis. Cell Death Dis. 11, 833. doi:10.1038/s41419-020-03038-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Ye, G., Xie, Z., Zeng, H., Wang, P., Li, J., Zheng, G., et al. (2020). Oxidative stress-mediated mitochondrial dysfunction facilitates mesenchymal stem cell senescence in ankylosing spondylitis. Cell Death Dis. 11, 775. doi:10.1038/s41419-020-02993-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, Y., Zhang, X., Chen, R., Han, S., Liu, Y., Liu, X., et al. (2020). Serum amino acid metabolic profiles of ankylosing spondylitis by targeted metabolomics analysis. Clin. Rheumatol. 39, 2325–2336. doi:10.1007/s10067-020-04974-z

PubMed Abstract | CrossRef Full Text | Google Scholar