lnc-SAMD14-4 can regulate expression of the COL1A1 and COL1A2 in human chondrocytes

Specimen collection

A Kellgren-Lawrence (KL) system was used to measure OA severity in patients (Kellgren & Lawrence, 1957). OA cartilage (KL grades 3 or 4) was collected from 34 patients undergoing total knee joint replacement due to severe OA. Normal cartilage (KL grades 0) was collected from 19 patients without history of rheumatoid arthritis or OA (who were being treated for trauma, thromboangiities obliterans, osteosarcoma, or limb amputation). Informed consent was signed by all patients and additional informed consent was obtained for publication of subject age and sex. This study was approved by the Ethics Commission of the 163rd Central Hospital of the Chinese People’s Liberation Army (IRB. [2015] 001). Patient data appears in Table S1. Clinical specimens were divided for qPCR and for array analyses.

LncRNA and mRNA microarray analysis

Quality testing, sample labeling, and hybridization to human WT lncRNA microarray (Affymetrix) were conducted by Oebiotech Co (Shanghai, China) (Fig. S1); their detailed protocol can be found in previous publications (Lightfoot, 2016; Mueller et al., 2015). All RNA integrity number (RIN) values were >7, and 260:280 ratios were >1.7 (as measured by NanoDrop). The Genechip included 63,542 non-coding and 27,134 coding transcripts. Images were taken with an Affymetrix Scanner3000, and data normalization and probe filtering were conducted using an Affymetrix GeneChip Expression Console (version 4.0; Affymetrix). The lncRNA-mRNA expression between OA and non-OA groups was identified as a ≥2.0 fold-change (FC) for up-regulation or a ≤ − 2.0 fold-change for down-regulation, with p ≤ 0.05. Afterwards, GO and KEGG analyses were used to identify roles of differentially expressed mRNAs with respect to GO terms or pathways. Their detailed definitions can be found in previous publications (Huang et al., 2011; Tao et al., 2012).

A volcano plot and hierarchical clustering were created to show distinctions in gene expression patterns among samples. Finally, a coding-non-coding gene co-expression (CNC) network was drawn using Cytoscape software. The microarray data are available through the GEO database with the accession number GSE113825.

RNA extraction and quantitative PCR

A RNeasy Lipid Tissue Mini Kit (Qiagen, Hilden, Germany) was used to extract total RNA from tissues according to the manufacturer’s directions. A NanoDrop ND-1000 spectrophotometer was used to measure RNA quality and purity. Then SDS-PAGE was used to measure RNA integrity. Total RNA from samples was reverse transcribed into cDNA using an All-in-one First-Strand cDNA Synthesis Kit (GeneCopoeia Inc, Rockville, MD, USA) as directed by the user manual. qPCR primers were designed and synthesized at the Shanghai Sangon Company (Table S2); qPCR was conducted using AceQ qPCR SYBR Green Master Mix (Vazyme Biotech Co, Jiangsu Sheng, China) on a Bio-Rad CFX ConnectTM Real-Time PCR detection system (Bio-Rad, Munchen, Germany) using 10 µl SYBR Master Mix, 0.4 µl forward primer (10 µM), 0.4 µl reverse primer (10 µM), 0.4 µl ROX Reference Dye 1, two µl cDNA, and 6.8 µl MiliQ water. The qPCR reaction conditions were as follows: pre-denaturation at 95 °C for 5 min, followed by 95 °C for 10 s, and 60 °C for 30 s for 40 cycles, and a final extension at 95 °C for 15 s, 60 °C for 60 s, and 95 °C for 15 s. The 2^−ΔΔCt method was employed to assess gene expression, as described in previous studies (Livak & Schmittgen, 2001). Each experiment was conducted in triplicate.

Cell culture and treatment

Human primary chondrocytes were procured from CHI Scientific, Inc (China). The cells were incubated in DMEM (Gibco, Waltham, MA, USA) medium supplemented with 14% FBS (Gibco, USA) in 6-well culture plates (6 × 10⁵∕well). Cell cultures were cultivated at 37 °C in a 5% CO₂ incubator and the medium was changed every three days. All studies were conducted using primary chondrocytes from passages 1-4. Cells were treated either with or without 10 ng/ml human recombinant interleukin-1 beta (IL-1 β) (Sigma-Aldrich, St Louis, MO, Germany) for 12 h, 24 h and 48 h.

Cell transfection

The overexpression vector, lnc-SAMD14-4(pcDNA3.1⁺-lnc-SAMD14-4), and blank vector, pcDNA3.1(pcDNA3.1⁺- vector), were designed by and purchased from GeneChem Co., Ltd. (Shanghai, China). The vectors were transfected into human primary chondrocytes using Lipofectamine 3000 (Invitrogen, Carlsbad, CA, USA). Three targeting siRNA, including si-lnc-SAMD14-4 (1#, 2# and 3#), and a control siRNA were designed and purchased from Ribo Biotechnology Co. Ltd. (Guangdong, China). These siRNA (100nM/well) were transfected into human primary chondrocytes to assess gain or loss of function of lnc-SAMD14-4.

Statistical analysis

Differences in expression of lncRNA-mRNA between OA and non-OA groups were assessed with a Student’s t test. GO and KEGG pathways were assessed using Fisher’s exact test. Other data were analyzed using one-way ANOVA and multiple groups were compared using Bonferroni’s method. GraphPad Prism 6.0 (GraphPad Software San Diego, CA) was applied for data analysis and drawing. Data are shown as means ± SD, and p < 0.05 was considered statistically significant.

Analysis of differential expression among lncRNAs-mRNAs in OA

We found 2,042 lncRNAs differed significantly in expression levels between OA and non-OA tissue (Table S3). The most upregulated lncRNA was lnc-SAMD14-4 and the most downregulated lncRNA was lnc-MSMP-2. The 15 most dysregulated (either down or up) lncRNAs are listed in Table 1.

mRNA expression data revealed dysregulation of 2,011 mRNAs in OA cartilage tissue, as compared to non-OA cartilage tissue; it was found that 1,386 mRNA were upregulated and 625 mRNAs were downregulated (Table S4). Among these mRNAs, the most upregulated protein coding RNA was MMP-13, which is responsible for regulating articular cartilage degradation (Fosang et al., 1996). MYOC was the most strongly down-regulated. The top 15 upregulated and downregulated mRNAs appear in Table 2. Data on differentially expressed genes (including lncRNA and mRNA) appear in Figs. 1A–1C.

Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of dysregulated mRNAs

The most enriched GOs targeted by upregulated mRNAs were extracellular matrix organization (GO:0030198), extracellular exosome (GO:0070062), and collagen binding (GO:0005518) in biological process, cellular process, and molecular function (Figs. 2A–2C, Tables S5–S7). Protein stabilization (GO:0050821), cilium (GO:0005929), and heparan sulfate proteoglycan binding (GO:0043395) were the most enriched GO pathways associated with downregulated transcripts (Figs. 2D–2F, Tables S8–S10).

We identified 27 pathways with upregulated mRNAs and the most enriched pathway was cytochrome P450 (Table S11). Additionally, 94 pathways corresponded to downregulated mRNAs, according to KEGG pathway analysis, and the top enriched pathway was focal adhesion (Table S12). The top 20 pathways enriched by downregulated and upregulated transcripts appear in Figs 3A and 3B.

LncRNA-mRNA co-expression network

Co-expression networks are of biological interest because co-expressed genes are controlled by the same transcriptional regulatory program, are functionally related, or are members of the same pathway or protein complex (Weirauch, 2011; Stuart et al., 2003). It has been confirmed that apoptosis is involved in OA pathogenesis (Zhang et al., 2016b; Hwang & Kim, 2015; Andrew et al., 2007; Kim & Blanco, 2007), and that the PI3K-Akt signaling pathway is important for apoptosis (Franke et al., 2003). We analyzed the 74 dysregulated mRNAs found to be enriched in the PI3K-Akt signaling pathway (Table S13) and identified 1,974 correlations for lncRNA-mRNA interaction pairs that satisfy p < 0.01 and Pearson correlation coefficients>0.95 (Table S14). Then, the top 500 of 1,974 lncRNA-mRNA pairs were used to draw the network (Table S15). Of the top 500 lncRNA-mRNA pairs, which included 57 mRNAs and 363 lncRNAs, 362 pairs had a positive association; the rest were negatively associated. It appears that a single mRNA may pair with multiple lncRNAs while most lncRNAs may only pair with a single mRNA (Fig. 4). Some lncRNA-mRNA pairs included the most dysregulated mRNA and lncRNA. For example, COL1A1 (mRNA, FC = 175.43) and COL1A2 (mRNA, FC = 85.32) were strongly associated with lnc-SAMD14-4 (lncRNA, FC = 86.11). All lncRNA-mRNA co-expressed pairs had Pearson correlation coefficients>0.98

Verification of lncRNAs and mRNAs expression by qRT-PCR

To confirm microarray profiling analysis, four lncRNAs (lnc-SAMD14-4, lnc-MARCKS-7, lnc-MSMP-2 and lnc-NPVF-4) and four mRNAs (SERPINF1, COL1A1, MYOC and PAK3) were selected and measured using qRT-PCR in 53 clinical specimens (34 OA and 19 non-OA specimens). The results from the qRT-PCR agreed with chip profiling analysis (Fig. 5).

Prediction of novel lncRNA targeting

We studied the most dysregulated lncRNA, lnc-SAMD14-4, for its potential function in OA development. To analyze the gene coding-potential of lnc-SAMD14-4, the open-reading frames (ORFs) and codon-substitution frequency (CSF) were measured according to former research protocol (Liu et al., 2014; Qu et al., 2016). The ORF Finder (National Center for Biotechnology Information) failed to predict ORFs more than 300 nt. Furthermore, BLAST search of all the positive-sense strand ORFs, which are marked by “ +” in Fig. S2, yielded no highly homologous protein among short peptides. Alternatively, PhyloCSF analysis of lnc-SAMD14-4 was negative (Fig. S3). These results confirmed that lnc-SAMD14-4 does not encode any functional protein product.

COL1A1, COL1A2 and lnc-SAMD14-4 were highly expressed in the IL-1 β-treated human primary chondrocytes

According to our co-expression network analysis, COL1A1 and COL1A2 were strongly associated with lnc-SAMD14-4. Thus, we speculated that lnc-SAMD14-4 may contribute to OA by interacting with COL1A1 and COL1A2. Previous studies have used human primary chondrocytes treated with IL-1 β-treated to model the process of human osteoarthritis (Fu et al., 2015; Goldring et al., 1988). We determined expression levels of COL1A1, COL1A2 and lnc-SAMD14-4 by qPCR. We found that COL1A1 and COL1A2 were upregulated after human chondrocytes were treated with IL-1 β for 24 h and 48 h. Furthermore, upregulation of lnc-SAMD14-4 with IL-1 β treatment was time-dependent (Fig. 6). These results suggest a positive correlation between lnc-SAMD14-4 and COL1A1, COL1A2 expression induced by IL-1 β treatment. Hence, human primary chondrocytes were treated either with or without IL-1 β for 48 h in the follow-up experiments.

lnc-SAMD14-4 was involved in IL-1 β-induced COL1A1 and COL1A2 expression in human primary chondrocytes

To further explore the function of lnc-SAMD14-4 in eliciting IL-1 β-induced COL1A1 and COL1A2 expression, we performed gain/loss of function experiments. We confirmed that lnc-SAMD14-4 expression in human primary chondrocytes increased after pcDNA3.1⁺- lnc-SAMD14-4 transfection (Fig. S4). Additionally, as shown in Fig. 7, lnc-SAMD14-4 overexpression promoted expression of COL1A1 and COL1A2 in the human primary chondrocytes, regardless of treatment with IL-1 β (p < 0.05).

To explore the knockdown effects of lnc-SAMD14-4 on COL1A1 and COL1A2, we designed and transfected three independent si-lnc-SAMD14-4 (1#, 2# and 3#) in human primary chondrocytes. The strongest silencing effects were observed in si-lnc-SAMD14-4 (2#), which was used for silencing in subsequent experiments (Fig. S5). The knockdown of lnc-SAMD14-4 significantly decreased the expression of COL1A1 and COL1A2 in human primary chondrocytes, regardless of treatment with IL-1 β-treatment (Fig. 8).