Screening of differentially expressed miRNAs and mRNAs in osteoarthritis.
After processing the raw data of mRNA expression profile GSE55457, GSE82107 and miRNA expression profile GSE143514, they were analyzed by the screening criteria |log2Fold change| > 1 and p < 0.05. Finally, 626 differentially expressed mRNAs were screened from GSE55457 (313 up-regulated and 313 down-regulated) (Fig. 1A and supplementary Fig. 1), 1647 differentially expressed mRNAs were screened from GSE82107 (498 down-regulated and 1147 up-regulated) (Fig. 1A and supplementary Fig. 1), and the same method was screened from GSE143514 to 142 Two differentially expressed miRNAs (22 down-regulated and 120 up-regulated) (Fig. 1A and supplementary Fig. 1). The Venn diagram was used to compare the differentially expressed mRNAs of the two expression profiles of GSE55457 and GSE82107, and 99 identical mRNAs were obtained (Fig. 1B).
Identifying Potential Mrnas Related To The Occurrence Of Osteoarthritis
In addition, in order to further identify the potential mRNAs in osteoarthritis, we re-downloaded an expression profile GSE55235 and performed WGCNA. First, the data of GSE55235 was standardized, and then 334 differentially expressed mRNAs (supplementary Fig. 2) were screened from the expression profile of GSE55235 through the limma package. Subsequently, I collected the expression levels of these differentially expressed genes for WGCNA. In order to ensure that the co-expression network is unsigned, we choose a soft threshold β = 16 (v 1C). Next, the average linkage hierarchical clustering method is used to cluster genes (Fig. 1D). At the same time, we perform cluster analysis on the modules and get a total of 3 modules (Fig. 1E). Then, according to the correlation between the characteristic genes of each module and the sample type, it is found that the turquoise module has the greatest correlation with the sample type (Fig. 1F). These results suggest that 27 genes in the turquoise module may be closely related to the occurrence of osteoarthritis. Moreover, we compared the 27 hub genes in the turquoise module with 99 intersections differentially expressed mRNAs, and a total of 19 intersection key mRNAs were obtained (Fig. 1G and Table 1).
Table 1
19 intersection key mRNAs from turquoise module and 99 intersections differentially expressed mRNAs.
|
mRNAs
|
19 intersection key mRNAs
|
APOLD1, CD52, CXCL13, CXCL2, CXCL3, DUSP1, FKBP5, FOSB, IGJ, IGLL5, JUN, MMP1, ZB1, NFIL3, NKG7, NR4A2, PLCD3, POU2AF1, TNFRSF17.
|
turquoise module
|
BCL6, RHOB, NFIL3, KLF4, PLCD3, CEBPD, GADD45A, BTG1, GADD45B, EIF1, APOLD1, CD52, NKG7, CXCL13, CXCL3, CXCL2, DUSP1, MMP1, JUN, FKBP5, NR4A2, FOSB, IGJ, TNFRSF17, POU2AF1, MZB1, IGLL5.
|
Functional Enrichment Analysis Of Potential Mrnas
In order to learn more about the hub mRNAs, we did the GO and KEGG analysis by Metascape. GO function enrichment analysis showed that 19 hub mRNAs mainly involved in immune related process. Just like what have showed in Fig. 2A, For example, the humoral immune response, cell response to chemokines, leukocyte chemotaxis, granulocyte chemotaxis, granulocyte migration, cell response to lipopolysaccharide, B cell activation and leukocyte migration, antibacterial humoral response as well as the defense responses against bacteria and so on. KEGG pathway enrichment analysis showed that 19 hub mRNAs were mainly involved in the IL-17 signaling pathway, TNF signaling pathway, TGF-beta signaling pathway, and nod-like receptor Signaling pathway, Chemokine signaling pathway, signaling by Interleukins and PI3K/AKT signaling pathway (Fig. 2B).
Ppi Network And Mirna-mrna Co‑expressed Network Analysis Of Hub Genes
Furthermore, in order to reveal the interaction relationship between these hub genes and find out the key genes further, we used String and Cytoscape to construct a PPI network (Fig. 2C), which consists of 11 nodes and 24 edges. Next, we obtained the highest scoring subnetwork in PPI through the MCODE plug-in, which was composed of five mRNAs including CXCL2, JUN, CXCL3, MMP1 and PLCD3(Fig. 2D). These genes are the most important genes in the PPI network and may be the most relevant biomarkers in osteoarthritis, so we selected these genes for further analysis. moreover, we used DIANA Tools online database to predict targeted miRNAs for 5 key genes (CXCL2, JUN, CXCL3, MMP1 and PLCD3). Then, all the targeted miRNAs obtained were compared with the differentially expressed miRNAs screened, 24 common miRNAs of 5 key genes were obtained (Fig. 2E). Also, 39 miRNA-mRNA pairs were predicted. According to the predicted results, a mRNA-miRNA co-expression network consisting of 29 nodes and 39 edges was constructed using Cytoscape (Fig. 2F and Table 2).
Table 2
miRNAs and mRNAs in miRNA-mRNA co-expression network.
mRNA
|
miRNAs
|
PLCD3
|
hsa-miR-1226-5p, hsa-miR-18a-3p, hsa-miR-21-3p, hsa-miR-34a-5p
|
MMP1
|
hsa-miR-205-5p, hsa-miR-20a-5p, hsa-miR-21-3p, hsa-miR-210-3p, hsa-miR-34a-5p, hsa-miR-520f-3p, hsa-miR-520c-3p, hsa-miR-941
|
CXCL3
|
hsa-miR-132-3p, hsa-miR-203a-3p, hsa-miR-20a-5p, hsa-miR-210-3p, hsa-miR-34a-5p, hsa-miR-520f-3p, hsa-miR-532-3p
|
JUN
|
hsa-miR-20b-5p, hsa-let-7e-5p, hsa-miR-138-5p, hsa-miR-144-3p, hsa-miR-17-5p, hsa-miR-18a-5p, hsa-miR-203a-3p, hsa-miR-20a-5p, hsa-miR-21-3p, hsa-miR-34a-5p, hsa-miR-490-3p, hsa-miR-93-5p, hsa-miR-93-3p
|
CXCL2
|
hsa-miR-138-5p, hsa-miR-15b-3p, hsa-miR-20a-5p, hsa-miR-21-3p, hsa-miR-210-3p, hsa-miR-215-5p, hsa-miR-34a-5p
|
Validation Of Hub Genes And Immune Infiltration Analysis
For the 5 hub genes obtained, the relative expression was detected by qRT-PCR to verify the analysis results. As shown in Fig. 3A. The relative expression level of PLCD3 was significantly increased in 15 pairs of synovial samples of osteoarthritis (p < 0.05), while the relative expression levels of MMP1 was significantly decreased (p < 0.05), which was consistent with our expected results. However, the qPCR results of JUN, CXCL3 and CXCL2 were inconsistent with the analysis (Table 3), which may be due to the deviation caused by the small number of samples. Furthermore, based on the expression profile data of GSE55457 and GSE82107, we used GraphPad Prism to analyze the clinical diagnostic value of 2 specifically expressed hub mRNAs (MMP1 and PLCD3) in osteoarthritis samples and draw the ROC curve. Compared with control samples, these 2 specifically expressed hub genes have higher diagnostic value in osteoarthritis samples. The AUC values were 0.7994 (MMP1, p < 0.001, Fig. 3B) and 0.7607 (PLCD3, p < 0.001, Fig. 3C), respectively. Therefore, we speculate that MMP1 and PLCD3 may be biomarkers for diagnosis of osteoarthritis. Among the pathway enrichment, MMP1 was mainly involved in the IL-17 signaling pathway, rheumatoid arthritis and leukocyte migration. PLCD3 was mainly involved in PI3K/AKT signaling pathway (Fig. 2B). At present, there are relatively few studies on PLCD3 in osteoarthritis, so we choose PLCD3 and PI3K/AKT pathways for follow-up research.
To confirm the role of PLCD3 in osteoarthritis, we used GSE82107 expression profile samples as the validation group. First, due to the importance of immune cells in the development of osteoarthritis, I analyzed the immune infiltration between osteoarthritis and normal tissue using CIBERSORT. Analysis results show that the proportions of B cells naïve, NK cells activated, Macrophages M2, Mast cells activated in osteoarthritis tissues were relatively high, and statistically significant(p < 0.05) (Fig. 3D). we investigated the correlation between the expression level of PLCD3 and above 4 immune cells (Fig. 3E and Fig. 3F). The results showed that the proportion of 4 immune cells was positively correlated with the expression of PLCD3.
Table 3
The fold change of PLCD3, MMP1, CXCL3, JUN and CXCL2 in GSE55457, GSE82107 and GSE55235.
mRNAs
|
Fold change(log2)
|
GSE55457
|
GSE82107
|
GSE55235
|
PLCD3
|
3.0550
|
2.9338
|
1.8978
|
MMP1
|
-1.48205
|
-5.4800
|
-1.41045
|
CXCL3
|
2.2569
|
4.5952
|
1.057053
|
JUN
|
1.7633
|
1.9807
|
1.6182
|
CXCL2
|
2.0021
|
1.56015
|
1.1196
|
PLCD3 is the direct target of miR-34a-5p.
First, the targeted miRNA miR-34a-5p, which was most correlated with PLCD3, was screened from the miRNA-mRNA co-expression network, and the qPCR detection showed that miR-34a-5p was low expressed in the synovial tissue of osteoarthritis, which was consistent with the biochemistry analysis results (Fig. 4A). Correlation analysis showed that the expression levels of miR-34a-5p and PLCD3 were negatively correlated (Fig. 4B). Then, the relationship between PLCD3 and miR-34a-5p was verified by dual luciferase assay. The complementary region sequence between 3-UTR PLCD3 and miR-34a-5p was shown in Fig. 4C. After transfection with miR-34a-5p mimic, the relative luciferase activity of PLCD3 was significantly decreased in WT group, but not in PLCD3 mut group (Fig. 4D). Meanwhile, qPCR analysis (Fig. 4E and Fig. 4F) showed that overexpression of miR-34a-5p could inhibit the expression of PLCD3, and knockdown of miR-34a-5p could promote the expression of PLCD3. In addition, the promotion effect of pCDH-PLCD3 on the expression of PLCD3 and the inhibition effect induced by si-PLCD3 were weakened by transfection of miR-34a-5p mimic or inhibitor, respectively. All the results suggested that PLCD3 was a direct target of miR-34a-5p.
The effect of miR-34a-5p/PLCD3 axis on the proliferation, migration, apoptosis and inflammatory factor production of HFLS-OA cells.
Subsequently, CCK8 assay (Fig. 4G) showed that si-PLCD3 significantly inhibited the proliferation of HOA -FLS cells. In contrast, the proliferation of HFLS-OA cells was significantly increased after transfection with pCDH-PLCD3. Co-transfection of miR-34a-5p mimic and pCDH-PLCD3 could change the corresponding effects of individual transfection of miR-34a-5p mimic and pCDH-PLCD3. Wound healing assay (Fig. 4H and Fig. 4I) showed an opposite trend. Overexpression of PLCD3 reduced cell migration and knockout of PLCD3 significantly increased cell migration rate. Similarly, co-transfection of miR-34a-5p mimic and pCDH-PLCD3 could reverse the effect of pCDH-PLCD3 on the migration ability of HFLS-OA cells. In conclusion, miR-34a-5p/PLCD3 axis may affect the proliferation and migration of HFLS-OA cells. To further clarify the effect of miR-34a-5p/PLCD3 axis in HFLS-OA cells, we performed western blotting. Western results showed that in HFLS-OA cells, overexpression of miR-34a-5p significantly increased the protein levels of Casp3 and Casp9, while overexpression PLCD3 showed the opposite trend (Fig. 5A and Fig. 5B). These results suggest that miR-34a-5p/PLCD3 axis can stimulate the apoptosis of HFLS-OA cells.
Furthermore, in order to determine the effect of miR-34a-5p/PLCD3 on the secretion of inflammatory factors in HFLS-OA cells, HFLS- OA cells were pretreated for 24 h with IL-1β (10 ng/ mL), miR-34a-5p mimic and PLCD3 overexpressed plasmid. And then TNFα and IL-6 content were detected with Elisa. As shown in Fig. 5C and Fig. 5D, in the single treatment group, the IL-1β group and PLCD3 overexpression group significantly increased the expression of TNFα and IL-6 compared with the control group. The miR-34a-5p mimic group showed protective activity. In the combined treatment group, PLCD3 overexpression promoted the pro-inflammatory cytokine secretion effect of IL-1β. On the contrary, miR-34a-5p mimic weakened the promotion of IL-1β on inflammatory factors.
miR-34a-5p /PLCD3 axis regulates HFLS-OA cells through PI3K/AKT signaling pathway.
In KEGG enrichment analysis, we found that PLCD3 was an enrichment gene in PI3K/AKT pathway. Therefore, whether miR-34a-5p/PLCD3 axis could regulate HFLS-OA cells through PI3K/AKT pathway? Western blot analysis was performed. As shown in Fig. 5G and Fig. 5H, overexpression of miR-34a-5p reduced the protein level of p-PI3K and p-AKT, while overexpression of PLCD3 showed the opposite trend of p-PI3K and p-AKT protein expression. Furthermore, HFLS-OA cells were treated and analyzed using 6-Bromoindirubin-3'-oxime (BIO), which was an inhibitor of the PI3K/AKT pathway. Firstly, different concentrations(0 µM, 0.625 µM, 1,25 µM, 2.5 µM, 5 µM, 10 µM, 20 µM and 50 µM) of BIO were set for cell activity experiment, and the IC50 value (IC50 = 5.95 µM) of BIO in HFLS-OA cells was obtained by fitting IC50 curve (Fig. 5E and Fig. 5F), which was used for subsequent interference experiments. The results showed that overexpression of miR-34a-5p increased the inhibitory effect of BIO on p-PI3K and p-AKT protein expression, while overexpression of PLCD3 significantly reversed the effect of BIO on p-PI3K and p-AKT protein expression (Fig. 5G and Fig. 5H). Therefore, we speculate that miR-34a-5p /PLCD3 may inhibit the activity of PI3K/AKT signaling pathway in HFLS-OA cells.