Jiaming Zhang
Xiaoxia Hao
Ruimin Chi3- 1Department of Orthopedics, Tongji Hospital, Tongji Medical College, Huazhong University of Science and Technology, Wuhan, China
- 2Cancer Center, Union Hospital, Tongji Medical College, Huazhong University of Science and Technology, Wuhan, China
- 3Department of Rehabilitation, Tongji Hospital, Tongji Medical College, Huazhong University of Science and Technology, Wuhan, China
Background: Osteoarthritis (OA) is a common degenerative joint disease. The aims of this study are to explore the effects of mechanical stress on whole transcriptome landscape and to identify a non-coding transcriptome signature of mechanical stress.
Methods: Next-generation RNA sequencing (RNA-seq) was performed on IL-1β-induced OA-like chondrocytes stimulated by mechanical stress. Integrated bioinformatics analysis was performed and further verified by experimental validations.
Results: A total of 5,022 differentially expressed mRNAs (DEMs), 88 differentially expressed miRNAs (DEMIs), 1,259 differentially expressed lncRNAs (DELs), and 393 differentially expressed circRNAs (DECs) were identified as the transcriptome response to mechanical stress. The functional annotation of the DEMs revealed the effects of mechanical stress on chondrocyte biology, ranging from cell fate, metabolism, and motility to endocrine, immune response, and signaling transduction. Among the DELs, ∼92.6% were identified as the novel lncRNAs. According to the co-expressing DEMs potentially regulated by the responsive DELs, we found that these DELs were involved in the modification of immune and metabolism. Moreover, immune- and metabolism-relevant DELs exhibited a notable involvement in the competing endogenous RNA (ceRNA) regulation networks. Silencing lncRNA TCONS_00029778 attenuated cellular senescence induced by mechanical stress. Moreover, the expression of Cd80 was elevated by mechanical stress, which was rescued by silencing TCONS_00029778.
Conclusion: The transcriptome landscape of IL-1β-induced OA-like chondrocytes was remarkably remodeled by mechanical stress. This study identified an immune- and metabolism-related ncRNA transcriptome signature responsive to mechanical stress and provides an insight of ncRNAs into chondrocyte biology and OA.
Introduction
Osteoarthritis (OA) is a common degenerative joint disease characterized by degradation of articular cartilage and subchondral bone, which causes a significant socioeconomic burden (Felson, 2010; Glyn-Jones et al., 2015). Mechanical stimulation from daily activities and physical exercise is indispensable for articular cartilage fitness. Physical exercise also gives promise to OA management by alleviating symptoms and structural progress (Ramage et al., 2009; Iijima et al., 2016; Chang et al., 2017). However, excessive mechanical loading is a critical risk factor of OA resulting in the damaged integrity of cartilage-subchondral bone unit, characterized by chondrocyte death, cartilage degeneration, and subchondral bone remodeling (Sharma et al., 2013; Poulet et al., 2015; Jorgensen et al., 2017).
As a key component of cartilage, chondrocytes are mechanosensitive, perceiving and responding to mechanical stress throughout life in a magnitude-dependent manner (Liu et al., 2016b). It is still a puzzle how the mechanical stimulation favors the cartilage and chondrocytes and how this process is attenuated and the outcome is turned over by excessive loading. In principle, moderate mechanical stress favors chondrocytes by regulating multiple signaling pathways and cellular functions, such as cytoskeleton and primary cilia (Yang et al., 2016; Xiang et al., 2018; Fu et al., 2019), while excessive mechanical stimulation always triggers the unfavorable changes and outcomes of chondrocytes, for instance, enhanced catabolic effects, compromised respiratory function, and apoptosis (Sun, 2010; Coleman et al., 2016). Mitochondrial dysfunction is a well-accepted hallmark of OA (Buckwalter et al., 2013). Mitochondrial fitness is critical for chondrocyte fate and functions, at least, by regulating metabolism (Zhang et al., 2021). Under unfavorable stress, chondrocytes not only increase reactive oxygen species (ROS) production, but also undergo metabolism shift from oxidative phosphorylation to glycolysis via the AMP-activated protein kinase (AMPK) and mechanistic target of rapamycin (mTOR) pathways (Zheng et al., 2021). Moreover, mitochondrial mechanotransduction is also critical for chondrocyte biology in the pathogenesis of OA. Activating by ROS/Liver kinase B1 (LKB1) and Ca2+/calcium/calmodulin-dependent kinase kinase 2 (CAMKK2), AMPK signaling pathway and subsequent ECM production, cytoskeleton rearrangement, and mitochondrial fitness. AMPK signaling pathway may be at the core of the events in which mechanical stimulation elicits alterations in mitochondrial functions, while the mechanisms of reduced AMPK activity in OA are elusive (Jiang et al., 2021).
Due to the application of high-throughput methods, increasing numbers of disturbed signaling pathways and mechanosensitive genes are highlighted, including long-noncoding RNAs (lncRNA). lncRNAs are broadly defined as a class of non-coding transcripts more than 200 nucleotides in length without coding potentiality (Zhang J et al., 2018). Liu et al. (2016a) identified that 107 lncRNAs were differentially expressed in the damaged cartilage and increased lncRNA MRS was responsive to cyclic tensile strain and involved in extracellular matrix (ECM) degradation by disrupting cytoskeleton and increasing matrix metalloproteinases (MMPs). Qiu et al. (2021) identified 16 upregulated and two downregulated miRNAs, as well as 44 upregulated lncRNAs and 39 downregulated lncRNAs in rat femur and tibia under mechanical stress. These lines of evidence suggest the underlying involvement of lncRNAs in mechanotransduction. However, the comprehensive understanding of mechanical response of chondrocyte transcriptome landscape, including mRNAs, microRNAs (miRNAs), lncRNAs, and circular RNAs (circRNAs), is still absent. The roles of noncoding RNAs (ncRNAs) in mechanotransduction remain elusive, while some ncRNAs were identified as the indispensable participators in OA pathogenesis (Li et al., 2018; Malemud, 2018; Kong et al., 2021).
In this study, we employed next-generation sequencing (RNA-seq) to explore the effects of mechanical stress on whole transcriptome landscape and to identify a non-coding transcriptome signature of mechanics response. This study identified the immune- and metabolism-related lncRNAs that may be involved in the mechanisms upon overloading and provides an insight of ncRNAs into chondrocyte biology and OA.
Materials and Methods
Isolation and Culture of Chondrocytes
Chondrocytes were isolated from the knee articular cartilage of 3-day-old Sprague–Dawley rats. The animal protocol was approved by the Experimental Animal Ethics Committee of Tongji Medical College, Huazhong University of Science and Technology (Wuhan, China). Briefly, cartilage tissue was minced and digested with 0.25% trypsin (Gibco, United States) for 30 min and 0.25% collagenase II (Invitrogen, United States) at 37°C for 6 h sequentially. Chondrocytes were collected and cultured in an incubator containing 5% CO2 at 37°C with Dulbecco’s DMEM/F12 (HyClone, United States) supplemented with10% fetal bovine serum (FBS, Gibco, United States) and antibiotics (100 U/ml penicillin G sodium, 100 μg/ml streptomycin sulfate) (Gibco, United States).
IL-1β Inducement on Rat Chondrocyte
It is well-documented that IL-1β can induce OA-like phenotypes in cultured chondrocytes, characterized by the reduction of anabolism (decreased ECM components) and the increase of catabolism (increased metalloproteinases) (Yao et al., 2019; Yao et al., 2021; Zhang et al., 2021). The concentrations of IL-1β ranging from 1 ng/ml to as much as 1,000 ng/ml, compared to <10 pg/ml in body fluids, are considered as supra-physiologic concentrations. Therefore, supra-physiologic concentrations of IL-1β (5–10 ng/ml) are required to induce OA-like changes in cartilage explants and chondrocyte cultures. In our study, chondrocytes were induced by recombinant rat IL-1β (5 ng/ml) (Catalog No. 400-01B, Peprotech, United States) for 24 h to obtain OA-like phenotypes (Yao et al., 2019).
Mechanical Stress Stimulation
Mechanical stress was applied to IL-1β-induced chondrocytes by using the four-point bending system as our previous studies described (Xu et al., 2013). Briefly, 1×106 chondrocytes at passage first were placed on a special vinyl plate (7.5 cm × 4 cm) and cultured until reaching 80% confluence. Then, the vinyl plate was transferred to a bending dish filled with 20 ml of DMEM/F12 medium in an incubator containing 5% CO2 at 37°C. Chondrocytes were exposed to mechanical stress with intensity of 2,000 μstrain at a frequency of 1 Hz for 4 h. For the control group, cells were cultured in the same condition without mechanical stress stimulation for 4 h.
RNA Isolation and Real-Time Polymerase Chain Reaction (RT-PCR)
Total RNA was extracted immediately after IL-1β intervention (24 h) and mechanical stress stimulation (4 h). Total RNAs were extracted using RNeasy kit (Qiagen, Germany) according to the manufacturer’s instructions. Reverse transcription (RT) was performed using QuantiTect Reverse Transcription kit (Qiagen, Germany). RT-PCR was performed using CFX Connect Real-Time System (Bio-Rad, Singapore) with SYBR green master mix (TOYOBO, Japan) as previously described (Yao et al., 2019). The thermocycling conditions were as follows: 30 s of polymerase activation at 95°C, followed by 40 cycles of 95°C for 5 s and 60°C for 30 s. Gene expression was normalized to Gapdh and calculated according to the 2−ΔΔCT method. The primers used in this study are listed in Supplementary Table S1. Each experiment was repeated three times.
Small Interfering RNA Silencing
The small interfering RNA (siRNA) used in this study was purchased from Ribbio Biological (Guangzhou, China). siRNA-TCONS_00029778: GCACCAGAAATGCCATGAA. Transfection was preformed using Lipofectamine 3000 transfection reagent according to manufacturer’s protocol (Guo et al., 2020). Forty-eight hours after transfection, the cells were subjected to IL-1β inducement and mechanical stimulation.
β-Galactosidase Staining
β-Galactosidase staining was performed to assess cellular senescence of the mechanical stress-stimulated chondrocytes. The staining was performed according to the manufacturer’s instructions (Senescence β-Galactosidase Staining Kit #9860, Cell Signaling, United States) (Itahana et al., 2007). The average number of senescent chondrocytes quantified in four independent fields selected by a blind way was used to represent one individual experiment. Five individual experiments were repeated.
Mitochondrion and Lysosome Staining
The subcellular location and morphology of mitochondrion and lysosome was accessed by Mito-Tracker Green (Beyotime, China) and Lyso-Tracker Red (Beyotime, China) staining according to the manufacturer’s instructions (Lin et al., 2021; Zhang et al., 2021). Lysosome number of 50 chondrocytes were quantified in a double-blind way by two individual researchers. The mean values were calculated by the data of 50 cells in an individual experiment. Five individual experiments were repeated.
Next Generation RNA-Sequencing (RNA-Seq)
Total RNA was isolated by Trizol (Invitrogen, United States) according to the manufacturer’s instructions. Each group contains three biological replicates. Library construction followed the standard protocols (Djebali et al., 2012) and was sequenced on the Illumina HiSeq 4000 platform by Genedenovo Biotechnology Co., Ltd (Guangzhou, China). We used fastp tool (v0.18.0) (Chen et al., 2018) to filter for high-quality reads by trimming adapters and removing reads that contained >10% unknown nucleotides, and reads with >50% low quality (Q-value ≤ 20) bases. Clean paired-end reads were then mapped to the Rattus norvegicus reference genome (Ensembl release 91 Rnor 6.0) with HISAT2 v.2.1.0 (Kim et al., 2015) with default parameters. The resulting files were sorted using SAMtools (v1.1.1) and subjected to HTSeq (0.9.0) to obtain the counts of each gene.
Raw count was filtered by a minimum expression threshold of 5 for coding genes and 2 for non-coding genes. The genes with low count (<5 or <2) in both groups were removed. The principal component analysis (PCA) was performed to determine the overall difference between the mechanical stress-stimulated group and control (David and Jacobs, 2014). Differential expression analysis of gene abundances was performed by R package edgeR (v3.11) (Nikolayeva and Robinson, 2014). Differentially expressed transcripts were defined by the criteria of |log2 (fold change)| > 1 (|log2 FC| > 1) and FDR < 0.05 for mRNA and lncRNA, p-value < 0.05 for miRNA and circRNA, respectively. The expression of transcripts was visualized by the R package pheatmap (v1.0.12) (Li et al., 2020). Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed by the R package clusterProfiler (v3.11) (Yu et al., 2012). Also, Gene Set Enrichment Analysis (GSEA) was employed to determine whether an a priori defined set of genes shows statistically significant differences between the groups or samples by GSEA software (v4.1.0) (Subramanian et al., 2005).
Novel lncRNA Identification
Th software Coding-Non-Coding Index (CNCI) (v2) (Sun et al., 2013) and Coding Potential Calculator (CPC) (Kong et al., 2007) were used to assess the protein-coding potential of novel transcripts by default parameters. Meanwhile, novel transcripts were mapped to SwissProt (http://www.expasy.ch/sprot/) sequence database to provide a high level of protein annotation (Bairoch and Apweiler, 2000). The intersection of both non-protein-coding potential results and non-protein annotation results were defined as the novel lncRNAs.
Regulation Mechanism Analysis and Function Annotation of lncRNAs
Antisense lncRNAs potentially regulate the host gene expression by silencing transcription and mRNA stability. In order to reveal the interaction between antisense lncRNA and mRNA, the software RNAplex (v2.4.16) was used to predict the complementary correlation of antisense lncRNA and mRNA based on the calculation of minimum free energy through thermodynamics structure (Tafer and Hofacker, 2008). Also, lncRNAs regulate their neighboring genes on the same allele in a cis-regulation mechanism. lncRNAs in less than 100 kb up-/downstream of a gene are likely to be cis-regulators. Moreover, lncRNAs regulate the expression of non-adjacent genes in a trans-regulation mechanism (Kopp and Mendell, 2018). The protein-coding genes that are potentially regulated by lncRNAs in antisense, cis-, or trans-mechanisms and also co-expression with corresponding lncRNAs were defined as the target genes of lncRNAs. The correlation between lncRNA and mRNA was calculated using the Pearson correlation coefficient and p-value. The Pearson correlation value was >0.9 or <−0.9, and p-value < 0.05 indicated that there could be a significant correlation between lncRNA and mRNA. The target genes of lncRNAs were subjected to enrichment analysis of GO functions and KEGG pathways.
CeRNA Network Mining and Function Annotation
Before ceRNA network construction, the targets of DEMIs were predicted by RNAhybrid (v2.1.2)/svm_light (v6.01), Miranda (v3.3a), and TargetScan (v7.0). Ritchie (2017) define the miRNA-mRNA/lncRNA/circRNA pairs based on complementary bases. The DEMs potentially regulated by DECs were defined based on expression correlation evaluated by Pearson correlation coefficient.
The principles of ceRNA network construction were strong binding potential and expression relevance. Expression correlation between DEMIs and DEMs/DELs/DECs was evaluated by the Pearson correlation coefficient (R). Pairs with R < −0.9 were selected as co-expressed miRNA-mRNA/lncRNA/circRNA pairs, among which the pairs with positively co-expressed mRNA-lncRNA/circRNA were defined as ceRNA pairs. To increase the accuracy of prediction, hypergeometric cumulative distribution function test was performed as previous study described (Xu et al., 2016). p-value less than 0.05 was considered significant. To functionally annotate the ceRNA network, GO enrichment and KEGG pathway analysis of the DEMs in the ceRNA network were performed. Transcript connectivity in the ceRNA network was defined as number of co-expressed targeted miRNAs, which was employed to assess the significance of genes in biological networks. The ceRNA network was visualized by the R package Circlize (Version 0.4.13) (Gu et al., 2014).
Statistical Analysis
The statistics methods used in this study include unpaired two-tailed Student’s t-test and one-way ANOVA with Dunnett’s multiple comparisons test. p < 0.05 was considered statistically significant. The application of each statistics method was specified in figure legends. GraphPad Prism 7.0 was used for the statistical analysis, and data were presented as means ± SD.
Results
Hallmark Remodeling of RNA Whole Transcriptome Landscape in IL-1β-Induced Rat OA-like Chondrocytes as a Response to Mechanical Stress
Given that the mechanical stimulation employed in this study induced the genes relevant to OA-like phenotypes, including Mmp9 and Mmp13 (encoding MMP9 and MMP13), we found that this stimulation promoted chondrocyte catabolic metabolism (Supplementary Figure S1A). Also, senescent chondrocytes were increased under this stimulation (Supplementary Figure S1B). RNA-seq analysis showed a distinct whole transcriptome profile of mechanical stress-stimulated OA-like chondrocytes, which included a total of 5,022 differentially expressed mRNAs (DEMs), 88 differentially expressed miRNAs (DEMIs), 1,258 differentially expressed lncRNAs (DELs), and 393 differentially expressed circRNAs (DECs) (Supplementary Table S2), and the extensive crosstalk among these four compositions by co-expressing relation (Figures 1A,B). Principal component analysis based on mRNA or lncRNA profiles showed the high relevance of lncRNA transcriptome with mRNA in IL-1β-induced rat OA-like chondrocytes (Figure 1C).
FIGURE 1. Hallmark remodeling of RNA whole transcriptomic landscape in OA-like chondrocytes as a response to mechanical stress. (A) Transcriptomic landscape of mechanical stress-stimulated OA-like chondrocytes, including mRNA, miRNA, and lncRNA/circRNA, based on gene co-expression analysis. (B) A total of 5,022 differentially expressed mRNAs (DEMs), 88 differentially expressed miRNAs (DEMIs), 1,259 differentially expressed lncRNAs (DELs), and 393 differentially expressed circRNAs (DECs) were identified in the comparison of the mechanical stress-stimulated group and the control group (FDR < 0.05, |log2FC| > 1). (C) Principal component analysis of the mechanical stress-stimulated group and the control group based on mRNA profiles or lncRNA profiles. The mRNA transcriptome shows high relevance with the lncRNA transcriptome of OA-like chondrocytes.
The mRNA Transcriptome Reveals the Critical Mechanical Stress-Responsive Genes and Signaling Pathways
To understand the response of IL-1β-induced rat OA-like chondrocytes to mechanical stress, we first focused on the mRNA transcriptome and defined the top 15 up- or downregulated genes in the rank of FDR or log2FC as the outstanding responsive genes (Figures 2A,B), among which Srebf1 (Sterol regulatory element-binding transcription factor 1), Trp53inp1 (Tumor protein P53 inducible nuclear protein 1), Slc30a2 (Solute carrier family 30 member 2), Cyp1a1 (Cytochrome P450 family subfamily A member 1), Zc3h12a (Zinc finger CCCH-type containing 12A), Aqp8 (Aquaporin 8), Nisch (Nischarin), Cstf2 (Cleavage stimulation factor subunit 2), and Wbp4 (WW domain binding protein 4) were identified as the significant genes based on the rank of FDR and FC (Table 1, Supplementary Figure S2).
FIGURE 2. The mRNA transcriptome reveals the critical mechanical stress-responsive genes and signaling pathways. (A) Top 15 up- or downregulated genes in the rank of FDR. (B) Top 15 up- or downregulated genes in the rank of log2FC. (C) KEGG pathway analysis of DEMs by ORA method. The relevant pathways are assigned into five groups, namely, cell growth and death, cell metabolism, cell motility and community, endocrine response, and signaling transduction. The color and size of dot represent FDR and GeneRatio of the corresponding KEGG pathway, respectively. FDR, false discovery rate, p-value estimated by one-way ANOVA and adjusted by the Benjamini–Hochberg method. GeneRatio indicates the gene number ratio in each KEGG pathway. (D) KEGG pathway analysis of total gene expression profile by the GSEA method. The relevant pathways are assigned into seven groups, namely, cell growth and death, cell metabolism, cell motility and community, genetic information processing, endocrine response, immune response, and signaling transduction. The color and size of dot represent NES and FDR of the corresponding KEGG pathway, respectively. NES, normalized enrichment score. A positive NES indicates that genes over-represented in the gene set are mostly upregulated in the comparison of mechanical stress versus control and a negative NES instead indicated the opposite. FDR, false discovery rate, p-value estimated by one-way ANOVA and adjusted by the Benjamini–Hochberg method. (E) RT-PCR validation of the essential genes contributing to relevant signaling pathway responsive to mechanical stress, including Rock2, Aldh3a1, Vamp8, Slc30a1, Cyp1a1, Npr1, and Srebf1. *p <0.05, **p <0.01, two-tailed Student’s t-test, mechanical stress versus control. (F) Lysosome (red) and mitochondrion (green) staining with or without mechanical stimulation. Total number of lysosomes and the colocation events with mitochondria is quantified in each chondrocyte. ns, no significance; **p <0.01, two-tailed Student’s t-test; mechanical stress versus control. nc, nucleus.
TABLE 1. Notable differentially expressed coding genes responsive to mechanical stress.
GO analysis of DEMs showed that mechanical stimulation is highly relevant to the biological process, such as response to stress, cell death, and biosynthetic and metabolic process (Supplementary Figure S3). Moreover, mechanical stimulation regulated enzyme binding, kinase activity, GTPase activity, and oxidoreductase activity (Supplementary Figure S3), suggesting a connection between mechanical stress and cellular metabolism. KEGG pathway analysis by the ORA method showed that the mechanical stress-relevant pathways were assigned into five categories, including cell growth and death, cell metabolism, cell motility and community, endocrine response, and signaling transduction, among which p53 signaling pathway, mitophagy, focal adhesion, thyroid hormone signaling pathway, and Rap1 signaling pathway were the most significant pathway in each category, respectively (Figure 2C). Given that a large amount of DEMs were involved in PI3K-Akt (357 DEMs, not shown) and MAPK signaling pathway (304 DEMs, not shown), these two pathways were considered as widely disturbed by mechanical stress.
Given that the ORA method only considers a group of pre-selected genes by a specific criterion (FDR < 0.05 and log2FC < −1 or >1), we also employed the GSEA method to explore the pathways relevant to global change of transcriptome. Different with the pathway categories identified by the ORA method, which were considered as the robust response of signaling pathway, three more categories, including genetic information processing, endocrine, and immune response, were further identified (Figure 2D; Supplementary Figure S4,S5). From the category perspective, cell metabolism, immune, and endocrine response were almost enhanced, while cell growth and death, cell motility and community, genetic information processing, and signaling transduction were almost suppressed. Specifically, the functions of ribosome and lysosome and the cellular metabolism, such as arachidonic acid, tyrosine, and fatty acid metabolism, were potentially activated, while ubiquitin-mediated proteolysis was assumed to be inhibited. Also, cell cycle and relevant pathways, such as DNA replication, were suppressed by mechanical stress, which is consistent with the item “cellular senescence” in the ORA results. Regarding the signaling pathways, mTOR, TGF-beta, and ErbB signaling pathways were likely to be suppressed, while PPAR signaling pathway was potentially activated.
The ten genes relevant to the signaling pathways responsive to mechanical stress (Figures 1A, 2B), namely, Gls, Mcoln1, Actn4, Rock2, Aldh3a1, Vamp8, Slc30a1, Cyp1a1, Npr1, and Srebf1, were selected for further verification by RT-PCR. Besides Gls, Mcoln1, and Actn4, seven out of ten were validated. There was increased expression of Aldh3a1, Vamp8, Slc30a1, Cyp1a1, Npr1, and Srebf1 and downregulation of Rock2 under mechanical stress (Figure 2E). Moreover, we found that the lysosomes co-located with mitochondria were increased while there was no difference in the total number of lysosomes between the stimulated and the control groups (Figure 2F).
Novel lncRNAs Emerged as a Marked Response to Mechanical Stress
Among the 6,649 lncRNA transcripts we detected, a total of 3,345 novel lncRNAs were identified by software CNCI, CPC, and the SwissProt database, accounting for more than 50% in total (Figure 3A). The density plot of the gene expression of these lncRNAs showed that mechanical stress increased the relatively high expressed lncRNAs and decreased the relatively low expressed lncRNAs, including both the known and the novel lncRNAs. The novel lncRNAs were the most significantly altered composition responsive to mechanical stress (Figure 3B). Moreover, as shown in Figure 3C, the number of novel DELs in the novel group was notably higher than the known lncRNAs. In this context, we considered that the novel lncRNAs may be the key response to mechanical stress. As shown in Figure 3D, the volcano plot of DELs showed the top 10 up- or downregulated known and novel DELs ranked by FDR, respectively. The top 10 known or novel DELs are summarized in Table 2. Among them, the known lncRNA AC127756.1 (top 1) and Pvt1 (top 3) and the novel lncRNA TCONS_00028770 (top 1) and TCONS_00062413 (top 2) were further validated by RT-PCR (Figure 3E). To functionally annotate these DELs, the coding genes that are potentially regulated by mechanical stress-responsive lncRNAs in trans-, cis-, and antisense manners were identified and further involved in KEGG pathway analysis. As shown in Figure 3F and Supplementary Figure S6, DELs were highly involved immune-relevant pathways, such as primary immunodeficiency, intestinal immune network for IgA production, autoimmune thyroid disease, allograft rejection, and graft-versus-host disease, and metabolism pathways, such as metabolism of xenobiotics by cytochrome P450 and tyrosine metabolism. Also, the DEM–DEL pairs involved in these immune and metabolism pathways were plotted in Figures 3E,G. TCONS_00016589-Blnk, TCONS_00077973-Itgb7, TCONS_00063105-Itgb7, TCONS_00029778-Cd80, TCONS_00015850-Bach1, TCONS_00035994-Mt2A, TCONS_00072487-Nktr, TCONS_00032987-Arap3, TCONS_00036854-Galnt2, and TCONS_00010901-Ccl6 were identified as the immune-relevant lncRNA-coding gene pairs according to the correlation. Interestingly, silencing lncRNA TCONS_00029778 attenuated cellular senescence induced by mechanical stress (Supplementary Figure S7). Moreover, the expression of Cd80 was elevated by mechanical stress, which was recused by silencing lncRNA TCONS_00029778 (Figure 3H).
FIGURE 3. Novel lncRNAs emerged as a key response to mechanical stress. (A) A total of 3,345 novel lncRNAs are identified by software CNCI, CPC, and the SwissProt database. (B) The density plot of the gene expression of the known and the novel lncRNAs. Mechanical stress increases the relatively high expressed lncRNAs and decreases the relatively low expressed lncRNAs, including the known and the novel lncRNAs. The novel lncRNAs are the most significantly altered composition responsive to mechanical stress. (C) The numbers of the non-DELs and DELs in the known and the novel groups. The number of DELs in the novel group is notably higher than the known group. (D) The volcano plot of DELs with FDR < 0.05 and log2FC > 1 or < −1 in the comparison of mechanical stress versus control. The red bubble represents the known lncRNA. The blue bubble represents the novel lncRNA. The upregulated known and novel lncRNAs are labeled by bright red and dark red, respectively. The downregulated known and novel lncRNAs are labeled by bright blue and dark blue, respectively. (E) RT-PCR validation of the known lncRNA Pvt1 and AC127756.1 and the novel lncRNA TCONS_00028770 and TCONS_00062413. *p <0.05, **p <0.01, two-tailed Student’s t-test, mechanical stress versus control. (F) The KEGG pathway analysis of the coding genes potentially regulated by mechanical stress-responsive lncRNAs in trans-, cis-, and antisense manners by the ORA method. The color and size of dot represent FDR and GeneRatio of the corresponding KEGG pathway, respectively. FDR, false discovery rate, p-value estimated by one-way ANOVA and adjusted by the Benjamini–Hochberg method. GeneRatio indicates the gene number ratio in each KEGG pathway. (G) Co-expression of lncRNA–mRNA pairs in trans-, cis-, or antisense regulation. Blue and yellow bands represent the cis- and antisense pairs, respectively. Red and blue label represent mRNA and lncRNA, respectively. T58664, the abbreviation of lncRNA TCONS_00058664. (H) RT-PCR validation of the lncRNA–mRNA TCONS_00029778-Cd80. MS, mechanical stress. *p <0.05, two-tailed Student’s t-test; #p <0.05, ##p <0.01, one-way ANOVA.
TABLE 2. Notable differentially expressed long non-coding RNA responsive to mechanical stress.
Notable Involvements of lncRNAs in ceRNA Regulation Networks
Besides the mechanisms of cis-, trans-, and antisense regulation, lncRNAs are also involved in ceRNA regulation networks and regulate target coding genes by the crosstalk with miRNAs and circRNAs. In this context, we also constructed the potential ceRNA regulation networks based on DEMs, DELs, DEMIs, and DECs. In Figures 4A,B, the top 5 up- or downregulated DEMIs and top 15 up- or downregulated DECs were identified as the outstanding responsive DEMIs or DECs. Two miRNA-centric ceRNA regulation networks were constructed, including networks 1 and 2. In network 1, mRNAs, lncRNAs, circRNAs, and miRNAs accounted for 47.5%, 44.4%, 5.6%, and 2.6%, respectively. A decreased proportion of lncRNAs (∼38.96%) and an increased proportion of circRNAs (∼44.64%) happened in network 2, while the numbers of miRNAs of both networks were comparable (Figure 4C). Given that the proportions of lncRNAs in both ceRNA networks were largest among ncRNAs, we speculated that lncRNAs may play a more critical role and serve as the hubs in network. Moreover, we calculated the degree of each composition, defined as the number of nodes, in ceRNA networks (Figure 4D). Indeed, the numbers of lncRNA hubs with the high degree > or = 5 were higher than both mRNAs and circRNAs. Moreover, two miRNA clusters miR-322/miR-20a/miR-17/miR-21/miR-298 and miR-140/miR-216 were identified in the subnetworks (Figure 4E). Interestingly, 705 overlapped DEMs connected to DELs or DECs were defined as the commonly regulated genes by lncRNAs and circRNAs, accounting for >99.9% of the DEMs regulated by circRNAs, while 467 DEMs were defined as the lncRNA specifically regulated genes. To know the underlying functions of lncRNAs in ceRNA networks, these two groups of genes were subjected to KEGG pathway analysis (Figures 4G,H). The lncRNA specifically regulated genes were mainly involved in metabolism pathways, such as glycine, serine, and threonine metabolism (Supplementary Figure S8) and type II diabetes mellitus, pyrimidine metabolism, and selenocompound metabolism, while both groups were involved in the pathways relevant to stress, apoptosis, and cell cycle.
FIGURE 4. Notable involvement of lncRNAs in ceRNA regulation networks. (A) Top 5 up- or downregulated miRNAs (DEMIs) in the rank of p-value. (B) Top 15 up- or downregulated circRNAs (DECs) in the rank of FDR. (C) miRNA-centric ceRNA regulation networks constructed by DEMs (red), DEMIs (blue), DELs (green), and DECs (purple). In the ceRNA network 1, DEMs, DELs, DECs, and DEMIs account for 47.5%, 44.4%, 5.6%, and 2.6%, respectively. In the ceRNA network 2, DEMs, DELs, DECs, and DEMIs account for 59.5%, 27.1%, 8.1%, and 5.3%. (D) The degree distribution of DEMs, DELs, and DECs in ceRNA networks. The degree is defined as the number of nodes connected. The degree > or = 5 is defined as high connectivity in networks. The DELs with high connectivity are more than the DEMs, while the DELs with the degree lower than 5 are less than DEMs. (E) The essential miRNA-centric subnetworks. (F) The overlap of the DEMs connected to DELs or DECs. A total of 705 DEMs are defined as the commonly regulated genes by lncRNAs and circRNAs, accounting for >99.9% in the DEMs regulated by circRNAs. A total of 467 DEMs are defined as the lncRNA specifically regulated genes. (G) The KEGG pathway analysis of the lncRNA specifically regulated genes by the ORA method. The color and size of dot represent FDR and count of the genes assigned into the corresponding KEGG pathway, respectively. FDR, false discovery rate, p-value estimated by one-way ANOVA and adjusted by the Benjamini–Hochberg method. GeneRatio indicates the gene number ratio in each KEGG pathway. (H) The KEGG pathway analysis of the commonly regulated genes by the ORA method.
Discussion
To understand the underlying ncRNA mechanisms of mechanical stress in OA pathogenesis, we employed RNA-seq to define the ncRNA transcriptome remodeling, characterized by thousands of differentially expressed mRNAs, lncRNA, miRNAs, and circRNAs and their complicated regulation relations, which provides a full landscape of transcriptome response. Moreover, these responsive ncRNAs are potentially biomarkers and therapeutic targets that could be modified in OA therapy as well as cartilage regeneration, such as implantation of pre-modified chondrocyte based on these ncRNAs.
In our study, we first defined the differentially expressed transcripts in each composition of whole transcriptome. In spite of the fact that mRNA composition has the largest number of differentially expressed transcripts, the ratio of DELs in lncRNA composition is up to 18.9%, the largest compared to others (mRNA: 12.8%, miRNA: 5.8%, circRNA: 5.4%), suggesting that lncRNAs are more responsive to mechanical stress. We speculated more notable roles of lncRNAs in whole transcriptome compared to miRNAs or circRNAs. Due to this speculation, we focused on the lncRNA composition and the relationship with the others in data analysis. Interestingly, we found more than 3,000 lncRNAs without the known annotation, which we defined as the novel lncRNAs. Moreover, mechanical stress gives promise to the novel lncRNA expression, and the number of the novel DELs was more than 10-fold of the known DELs, which suggests that this part is more responsive to mechanical stress in the lncRNA transcriptome. Indeed, the roles of lncRNAs are highlighted in chondrocyte biology and OA pathogenesis (Tu et al., 2020; van Hoolwerff et al., 2020), such as lncRNA Pvt1 (Li et al., 2017). In spite of the fact that increased Pvt1 is proven to be a promoter of chondrocyte apoptosis contributing OA, we found that Pvt1 was decreased instead in our system, suggesting that apoptosis might not be robustly induced by this mechanical stress. Indeed, the apoptosis-relevant pathways were not significantly responsive as the others reported (Kong et al., 2013; Xu et al., 2019), while a group of apoptosis-relevant DEMs were involved in this process. Except for a unique form of OA resulting from clear trauma, namely, post-traumatic osteoarthritis (POTA), chondrocyte apoptosis is a final outcome in late OA, while the effect of mechanical stress lasts throughout the lifespan, which may not be a pivotal inducing factor but an accumulated influence on chondrocyte fate. At the same time, we also found decreased proliferation and blocked cell cycle, indicating that the stimulation in our system induced cell senescence. Indeed, the senescent chondrocytes were increased under this mechanic system.
Non-random positioning of lncRNAs in genome implies a connection between lncRNA function and the genomic neighborhood (Statello et al., 2021). lncRNAs regulate neighborhood in cis-acting regulation (locally), a special form of which depends on the interaction between antisense lncRNAs and host genes, namely, antisense-regulation (Villegas and Zaphiropoulos, 2015). Besides this cis-acting regulation, lncRNAs also regulate distal protein-coding genes in trans (non-locally) (Yan et al., 2017). The knowledge of lncRNA-mediated cis-, trans-, and antisense-regulation is critical for our understanding of RNA-mediated gene regulation and genome complexity of mechanical response. Before that, we first need to understand the response of coding genes and then build the connection of coding genes and lncRNAs to define the functions of mechanical-responsive lncRNAs. In the functional annotation of coding genes, we found that the robust responses of signaling pathways lay on enhanced cell metabolism and immune response. Up to now, no evidence supports the involvement of lncRNA in metabolism modification in mechanical stress-induced chondrocytes, while it is widely accepted that mechanical stimulation alters cartilage and chondrocyte metabolism (Zheng et al., 2021; Bleuel et al., 2015). Our study revealed some metabolism-relevant lncRNA-coding gene pairs, such as TCONS_00015386-Adcy5. Adcy5 encodes a member of the membrane-bound adenylyl cyclase enzymes involved in various signaling pathways, such as cGMP-PKG (Gorbe et al., 2010) and cAMP (Vatner et al., 2015) signaling pathways, suggesting that mechanics-responsive TCONS_00015386 may also involve these pathways.
Pereira et al. (2016) showed that chondrocytes involve immune response by directly affecting T-cell proliferation and indirectly inhibiting monocyte differentiation, indicating the native functions of chondrocyte in immune regulation and also the underlying immune-relevant mechanisms of chondrocytes in OA pathogenesis. Indeed, a strong relevance of the mechanics-responsive lncRNAs to metabolism and immune pathway suggests a potential involvement of lncRNAs in the response to mechanical stress. Immune-related lncRNAs are being identified during the activation of the innate immune response and T-cell development, differentiation, and activation (Heward and Lindsay, 2014). In our study, some immune-relevant lncRNAs-coding gene pairs were identified, such as TCONS_00029778-Cd80. Given that CD80 plays a critical role in the initiation of T-cell responses and Pereira et al. (2016) built the link between T cells and chondrocytes, we speculated that lncRNA-mediated CD80 elevation in mechanical-induced chondrocytes may be involved in the immune mechanisms of OA pathogenesis. Indeed, the development of antigen-induced arthritis was significantly suppressed by inhibition of CD80 (Odobasic et al., 2008). Also, Ji et al. (2019) identified a subpopulation of chondrocytes by single-cell sequencing, namely, regulatory chondrocytes (RegCs), which highly expresses CD74, CD80, CD86, and HLADPA1 and involves immune-relevant signaling pathway, such as Toll-like receptor, JAK/STAT, and chemokine signaling. This funding gives a novel layer of understanding that chondrocytes are directly involved in immune response. These lines of evidence provide new insights and potential therapeutic directions of targeting-CD80 and CD80-relevant lncRNAs in overloading-induced cartilage degeneration.
Besides the transcriptional regulations, lncRNAs regulate gene expression in ceRNA mechanisms at the post-transcriptional level (Dykes and Emanueli, 2017; He et al., 2019). Accumulating numbers of studies have demonstrated the ceRNA mechanism in chondrocyte biology and OA development (Liu et al., 2016c). In our study, the significant interaction of miRNAs, lncRNAs, and circRNA was identified as a complicated ceRNA regulatory network, in which lncRNAs showed more outstanding involvements compared to circRNAs. Indeed, some lncRNAs were proven to be critical regulators in ceRNA mechanism, such as lncRNA ADAMTS9-AS2 (Huang et al., 2019), DANCR (Zhang L et al., 2018), and FOXD2-AS1 (Cao et al., 2018). We found that lncRNAs specifically regulate amino acid (AA) metabolism pathways, such as glycine, serine, and threonine metabolism and pyrimidine metabolism, in ceRNA mechanisms. Alterations of AA metabolism, such as glutamate- and arginine-family AA as well as their related metabolites, have also been identified in OA (Li et al., 2016). At the same time, the involvements of miRNAs and circRNA are also indispensable in the ceRNA network. Among miRNA clusters we identified, miR-20a (Zhao and Gong, 2019), miR-21 (Sekar, 2020), miR-322 (Georgieva et al., 2020), and miR-17 (Hu et al., 2018) have been validated in OA pathogenesis and miR-298 is responsive to IL-1β in OA chondrocyte (Akhtar et al., 2010). Moreover, miR-140 has a critical role in cartilage development and homeostasis (Araldi and Schipani, 2010) as well as OA pathogenesis (Miyaki et al., 2009) and miR-216a-5p regulates the proliferation and apoptosis of chondrocytes in OA by targeting the JAK-STAT signaling pathway (Zhang L et al., 2018). These lines of evidence suggest the significance of these miRNA clusters and the underlying mechanisms of mechanical stress by regulating these miRNAs.
In our previous comparison study (Zhang et al., 2021), we found the different or same effects between moderate (2,000 μstrain) and extensive (5,000 μstrain) mechanical stress on chondrocyte mitochondrial functions. We observed the effects of mechanical stress on mitochondrial functions and cell fate, or more specifically, apoptosis. We found a beneficial effect of 2,000 μstrain on cell survival, while we did not consider the impact on anabolic and catabolic metabolisms of chondrocytes. The previous findings seem to be contradicted to the present study. Actually, cellular senescence is defined as an apoptosis-resistant state with increased senescence-associated secretory phenotype (SASP). Of note, Mmp13 is a well-characterized component of SASP. In the present study, we found increased senescent chondrocytes under 2,000 μstrain, suggesting more apoptosis-resistant cells as we found in the previous study. Moreover, we found that the co-location of mitochondria and lysosomes was increased under this stimulation, which suggests the potential enhanced mitophagy. This is coherent with our previous finding that 2,000 μstrain is able to promote mitochondrial dynamic and mitophagy. Taking these two studies together, we concluded that 2,000 μstrain mechanical stress can induce OA phenotypes partially in an in vitro model, characterized by increased MMPs and accumulated cellular senescence. Chondrocyte senescence appears to be a compromise for survival in mechanical stress. The distinct effects of this stress may lie on the different responses of subpopulations among the primary chondrocytes, while our findings are based on the bulk effects of overall chondrocytes.
The limitations of this study should be discussed carefully. Firstly, we employed IL-1β-treated rat/rodent OA-like chondrocytes but not human OA chondrocytes isolated from patients. Also, our conclusion is based on the 5 ng/ml IL-1β treatment because we did not employ the multiple concentrations of IL-1β. Secondly, we employed only one mechanical model in this study, while chondrocytes are assumed to be subjected to multiple kinds of mechanical stimulation in vivo. Thirdly, we did not include the normal chondrocytes treated with mechanical stimulation but not IL-1β. As a result, we did not conclude the effects of mechanical stimulation in normal chondrocytes. Whether the effects of mechanical stimulation on OA-like chondrocyte transcriptome is IL-1β-dependent or not is still inconclusive. Finally, the low homology between human and rodent lncRNAs would hammer the significance of this study, while the functions of some lncRNAs are assumed to be resistant to evolutionary constraints at both RNA sequence and structural levels due to their critical roles (Karner et al., 2020).
Conclusion
In conclusion, mechanical stress reprograms whole transcriptome of IL-1β-induced OA-like chondrocytes and drives disruptions of signaling pathways. The mechanics-responsive lncRNAs may be critical in mechanical stress-induced metabolism and immune modifications. Further research into the functions of these genes and their interaction with lncRNAs may lead to a better understanding of their biogenesis and the mechanism underlying OA, and thus, they may eventually be used as diagnostic biomarkers and therapeutic agents for OA.
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 in the article/Supplementary Material.
Ethics statement
The animal protocol is approved by the Experimental Animal Ethics Committee.
Author Contributions
JZ and XH, the conception and design of the study; RC, JL, XS, and XD, data gathering; JZ and XH, data analyses and interpretations; JZ, JQ, and TX, drafting and final approval of the article.
Funding
This study was supported by the National Natural Science Foundation of China (Grant No. 82072556).
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.821508/full#supplementary-material
Abbreviations
AA, amino acid; Aqp8, Aquaporin 8; Acan, aggrecan; CPC, Coding Potential Calculator; CNCI, Coding-Non-Coding Index; CeRNA, endogenous RNA; CircRNAs, circular RNAs; Col2a1, collagen II alpha1; Cstf2, Cleavage stimulation factor subunit 2; Cyp1a1, Cytochrome P450 family subfamily A member 1; DECs, differentially expressed circRNAs; DELs, differentially expressed lncRNAs; DEMs, differentially expressed mRNAs; DEMIs, differentially expressed miRNAs; ECM, extracellular matrix; GO, gene ontology; GSEA, gene set enrichment analysis; LncRNAs, long-noncoding RNAs; KEGG, Kyoto Encyclopedia of Genes and Genomes; Nisch, Nischarin; MiRNAs, microRNAs; MMPs, metalloproteinases; NcRNAs, noncoding RNAs; OA, osteoarthritis; PCA, principal component analysis; RT-PCR, real-time polymerase chain reaction; RT, reverse transcription; Srebf1, Sterol regulatory element-binding transcription factor 1; Trp53inp1, Tumor protein P53 inducible nuclear protein 1; Slc30a2, Solute carrier family 30 member; Wbp4, WW domain binding protein 4; Zc3h12a, Zinc finger CCCH-type containing 12A.
References
Akhtar, N., Rasheed, Z., Ramamurthy, S., Anbazhagan, A. N., Voss, F. R., and Haqqi, T. M. (2010). MicroRNA-27b Regulates the Expression of Matrix Metalloproteinase 13 in Human Osteoarthritis Chondrocytes. Arthritis Rheum. 62, 1361–1371. doi:10.1002/art.27329
Araldi, E., and Schipani, E. (2010). MicroRNA-140 and the Silencing of Osteoarthritis. Genes Dev. 24, 1075–1080. doi:10.1101/gad.1939310
Bairoch, A., and Apweiler, R. (2000). The SWISS-PROT Protein Sequence Database and its Supplement TrEMBL in 2000. Nucleic Acids Res. 28, 45–48. doi:10.1093/nar/28.1.45
Bleuel, J., Zaucke, F., Brüggemann, G.-P., and Niehoff, A. (2015). Effects of Cyclic Tensile Strain on Chondrocyte Metabolism: a Systematic Review. PLoS One 10, e0119816. doi:10.1371/journal.pone.0119816
Buckwalter, J. A., Anderson, D. D., Brown, T. D., Tochigi, Y., and Martin, J. A. (2013). The Roles of Mechanical Stresses in the Pathogenesis of Osteoarthritis. Cartilage 4, 286–294. doi:10.1177/1947603513495889
Cao, L., Wang, Y., Wang, Q., and Huang, J. (2018). LncRNA FOXD2-AS1 Regulates Chondrocyte Proliferation in Osteoarthritis by Acting as a Sponge of miR-206 to Modulate CCND1 Expression. Biomed. Pharmacother. 106, 1220–1226. doi:10.1016/j.biopha.2018.07.048
Chang, N.-J., Lee, K.-W., Chu, C.-J., Shie, M.-Y., Chou, P.-H., Lin, C.-C., et al. (2017). A Preclinical Assessment of Early Continuous Passive Motion and Treadmill Therapeutic Exercises for Generating Chondroprotective Effects after Anterior Cruciate Ligament Rupture. Am. J. Sports Med. 45, 2284–2293. doi:10.1177/0363546517704847
Chen, S., Zhou, Y., Chen, Y., and Gu, J. (2018). Fastp: an Ultra-fast All-In-One FASTQ Preprocessor. Bioinformatics 34, i884–i890. doi:10.1093/bioinformatics/bty560
Coleman, M. C., Ramakrishnan, P. S., Brouillette, M. J., and Martin, J. A. (2016). Injurious Loading of Articular Cartilage Compromises Chondrocyte Respiratory Function. Arthritis Rheumatol. 68, 662–671. doi:10.1002/art.39460
David, C. C., and Jacobs, D. J. (2014). Principal Component Analysis: a Method for Determining the Essential Dynamics of Proteins. Methods Mol. Biol. 1084, 193–226. doi:10.1007/978-1-62703-658-0_11
Djebali, S., Davis, C. A., Merkel, A., Dobin, A., Lassmann, T., Mortazavi, A., et al. (2012). Landscape of Transcription in Human Cells. Nature 489, 101–108. doi:10.1038/nature11233
Dykes, I. M., and Emanueli, C. (2017). Transcriptional and Post-transcriptional Gene Regulation by Long Non-coding RNA. Genomics Proteomics Bioinformatics 15, 177–186. doi:10.1016/j.gpb.2016.12.005
Felson, D. T. (2010). Identifying Different Osteoarthritis Phenotypes through Epidemiology. Osteoarthr. Cartil. 18, 601–604. doi:10.1016/j.joca.2010.01.007
Fu, S., Thompson, C. L., Ali, A., Wang, W., Chapple, J. P., Mitchison, H. M., et al. (2019). Mechanical Loading Inhibits Cartilage Inflammatory Signalling via an HDAC6 and IFT-dependent Mechanism Regulating Primary Cilia Elongation. Osteoarthr. Cartil. 27, 1064–1074. doi:10.1016/j.joca.2019.03.003
Georgieva, V. S., Etich, J., Bluhm, B., Zhu, M., Frie, C., Wilson, R., et al. (2020). Ablation of the miRNA Cluster 24 Has Profound Effects on Extracellular Matrix Protein Abundance in Cartilage. Int. J. Mol. Sci. 21, 4112. doi:10.3390/ijms21114112
Glyn-Jones, S., Palmer, A. J. R., Agricola, R., Price, A. J., Vincent, T. L., Weinans, H., et al. (2015). Osteoarthritis. Lancet 386, 376–387. doi:10.1016/s0140-6736(14)60802-3
Gorbe, A., Giricz, Z., Szunyog, A., Csont, T., Burley, D. S., Baxter, G. F., et al. (2010). Role of cGMP-PKG Signaling in the protection of Neonatal Rat Cardiac Myocytes Subjected to Simulated Ischemia/reoxygenation. Basic Res. Cardiol. 105, 643–650. doi:10.1007/s00395-010-0097-0
Gu, Z., Gu, L., Eils, R., Schlesner, M., and Brors, B. (2014). Circlize Implements and Enhances Circular Visualization in R. Bioinformatics 30, 2811–2812. doi:10.1093/bioinformatics/btu393
Guo, J., Ren, R., Sun, K., Yao, X., Lin, J., Wang, G., et al. (2020). PERK Controls Bone Homeostasis through the Regulation of Osteoclast Differentiation and Function. Cell Death Dis. 11, 847. doi:10.1038/s41419-020-03046-z
He, R.-Z., Luo, D.-X., and Mo, Y.-Y. (2019). Emerging Roles of lncRNAs in the post-transcriptional Regulation in Cancer. Genes Dis. 6, 6–15. doi:10.1016/j.gendis.2019.01.003
Heward, J. A., and Lindsay, M. A. (2014). Long Non-coding RNAs in the Regulation of the Immune Response. Trends Immunol. 35, 408–419. doi:10.1016/j.it.2014.07.005
Hu, J., Wang, Z., Shan, Y., Pan, Y., Ma, J., and Jia, L. (2018). Long Non-coding RNA HOTAIR Promotes Osteoarthritis Progression via miR-17-5p/FUT2/β-Catenin axis. Cel Death Dis. 9, 711. doi:10.1038/s41419-018-0746-z
Huang, M.-J., Zhao, J.-Y., Xu, J.-J., Li, J., Zhuang, Y.-F., and Zhang, X.-L. (2019). lncRNA ADAMTS9-AS2 Controls Human Mesenchymal Stem Cell Chondrogenic Differentiation and Functions as a ceRNA. Mol. Ther. Nucleic Acids 18, 533–545. doi:10.1016/j.omtn.2019.08.027
Iijima, H., Aoyama, T., Ito, A., Tajino, J., Yamaguchi, S., Nagai, M., et al. (2016). Exercise Intervention Increases Expression of Bone Morphogenetic Proteins and Prevents the Progression of Cartilage-Subchondral Bone Lesions in a post-traumatic Rat Knee Model. Osteoarthr. Cartil. 24, 1092–1102. doi:10.1016/j.joca.2016.01.006
Itahana, K., Campisi, J., and Dimri, G. P. (2007). Methods to Detect Biomarkers of Cellular Senescence. Methods Mol. Biol. 371, 21–31. doi:10.1007/978-1-59745-361-5_3
Ji, Q., Zheng, Y., Zhang, G., Hu, Y., Fan, X., Hou, Y., et al. (2019). Single-cell RNA-Seq Analysis Reveals the Progression of Human Osteoarthritis. Ann. Rheum. Dis. 78, 100–110. doi:10.1136/annrheumdis-2017-212863
Jiang, W., Liu, H., Wan, R., Wu, Y., Shi, Z., and Huang, W. (2021). Mechanisms Linking Mitochondrial Mechanotransduction and Chondrocyte Biology in the Pathogenesis of Osteoarthritis. Ageing Res. Rev. 67, 101315. doi:10.1016/j.arr.2021.101315
Jørgensen, A. E. M., Kjær, M., and Heinemeier, K. M. (2017). The Effect of Aging and Mechanical Loading on the Metabolism of Articular Cartilage. J. Rheumatol. 44, 410–417. doi:10.3899/jrheum.160226
Karner, H., Webb, C.-H., Carmona, S., Liu, Y., Lin, B., Erhard, M., et al. (2020). Functional Conservation of LncRNA JPX Despite Sequence and Structural Divergence. J. Mol. Biol. 432, 283–300. doi:10.1016/j.jmb.2019.09.002
Kim, D., Langmead, B., and Salzberg, S. L. (2015). HISAT: a Fast Spliced Aligner with Low Memory Requirements. Nat. Methods 12, 357–360. doi:10.1038/nmeth.3317
Kong, L., Zhang, Y., Ye, Z.-Q., Liu, X.-Q., Zhao, S.-Q., Wei, L., et al. (2007). CPC: Assess the Protein-Coding Potential of Transcripts Using Sequence Features and Support Vector Machine. Nucleic Acids Res. 35, W345–W349. doi:10.1093/nar/gkm391
Kong, D., Zheng, T., Zhang, M., Wang, D., Du, S., Li, X., et al. (2013). Static Mechanical Stress Induces Apoptosis in Rat Endplate Chondrocytes through MAPK and Mitochondria-dependent Caspase Activation Signaling Pathways. PLoS One 8, e69403. doi:10.1371/journal.pone.0069403
Kong, H., Sun, M. L., Zhang, X. A., and Wang, X. Q. (2021). Crosstalk Among circRNA/lncRNA, miRNA, and mRNA in Osteoarthritis. Front. Cell Dev. Biol. 9, 774370.
Kopp, F., and Mendell, J. T. (2018). Functional Classification and Experimental Dissection of Long Noncoding RNAs. Cell 172, 393–407. doi:10.1016/j.cell.2018.01.011
Li, Y., Xiao, W., Luo, W., Zeng, C., Deng, Z., Ren, W., et al. (2016). Alterations of Amino Acid Metabolism in Osteoarthritis: its Implications for Nutrition and Health. Amino Acids 48, 907–914. doi:10.1007/s00726-015-2168-x
Li, Y., Li, S., Luo, Y., Liu, Y., and Yu, N. (2017). LncRNA PVT1 Regulates Chondrocyte Apoptosis in Osteoarthritis by Acting as a Sponge for miR-488-3p. DNA Cel Biol. 36, 571–580. doi:10.1089/dna.2017.3678
Li, H. Z., Lin, Z., Xu, X. H., Lin, N., and Lu, H. D. (2018). The Potential Roles of circRNAs in Osteoarthritis: a Coming Journey to Find a Treasure. Biosci. Rep. 38, BSR20180542. doi:10.1042/BSR20180542
Li, X., Yang, Y., Sun, G., Dai, W., Jie, X., Du, Y., et al. (2020). Promising Targets and Drugs in Rheumatoid Arthritis. Bone Jt. Res. 9, 501–514. doi:10.1302/2046-3758.98.bjr-2019-0301.r1
Lin, X., Wang, H., An, X., Zhang, J., Kuang, J., Hou, J., et al. (2021). Baeckein E Suppressed NLRP3 Inflammasome Activation through Inhibiting Both the Priming and Assembly Procedure: Implications for Gout Therapy. Phytomedicine 84, 153521. doi:10.1016/j.phymed.2021.153521
Liu, Q., Hu, X., Zhang, X., Dai, L., Duan, X., Zhou, C., et al. (2016a). The TMSB4 Pseudogene LncRNA Functions as a Competing Endogenous RNA to Promote Cartilage Degradation in Human Osteoarthritis. Mol. Ther. 24, 1726–1733. doi:10.1038/mt.2016.151
Liu, Q., Hu, X., Zhang, X., Duan, X., Yang, P., Zhao, F., et al. (2016b). Effects of Mechanical Stress on Chondrocyte Phenotype and Chondrocyte Extracellular Matrix Expression. Sci. Rep. 6, 37268. doi:10.1038/srep37268
Liu, Q., Zhang, X., Hu, X., Dai, L., Fu, X., Zhang, J., et al. (2016c). Circular RNA Related to the Chondrocyte ECM Regulates MMP13 Expression by Functioning as a MiR-136 'Sponge' in Human Cartilage Degradation. Sci. Rep. 6, 22572. doi:10.1038/srep22572
Miyaki, S., Nakasa, T., Otsuki, S., Grogan, S. P., Higashiyama, R., Inoue, A., et al. (2009). MicroRNA-140 Is Expressed in Differentiated Human Articular Chondrocytes and Modulates Interleukin-1 Responses. Arthritis Rheum. 60, 2723–2730. doi:10.1002/art.24745
Nikolayeva, O., and Robinson, M. D. (2014). edgeR for Differential RNA-Seq and ChIP-Seq Analysis: an Application to Stem Cell Biology. Methods Mol. Biol. 1150, 45–79. doi:10.1007/978-1-4939-0512-6_3
Odobasic, D., Leech, M. T., Xue, J. R., and Holdsworth, S. R. (2008). Distinctin Vivoroles of CD80 and CD86 in the Effector T-Cell Responses Inducing Antigen-Induced Arthritis. Immunology 124, 503–513. doi:10.1111/j.1365-2567.2007.02802.x
Pereira, R. C., Martinelli, D., Cancedda, R., Gentili, C., and Poggi, A. (2016). Human Articular Chondrocytes Regulate Immune Response by Affecting Directly T Cell Proliferation and Indirectly Inhibiting Monocyte Differentiation to Professional Antigen-Presenting Cells. Front. Immunol. 7, 415. doi:10.3389/fimmu.2016.00415
Poulet, B., De Souza, R., Kent, A. V., Saxon, L., Barker, O., Wilson, A., et al. (2015). Intermittent Applied Mechanical Loading Induces Subchondral Bone Thickening that May Be Intensified Locally by Contiguous Articular Cartilage Lesions. Osteoarthr. Cartil. 23, 940–948. doi:10.1016/j.joca.2015.01.012
Qiu, Y., Zhu, G., Zeng, C., Yuan, S., Qian, Y., Ye, Z., et al. (2021). Next Generation Sequencing of miRNAs and lncRNAs from Rat Femur and Tibia under Mechanical Stress. Mol. Med. Rep. 24, 561. doi:10.3892/mmr.2021.12200
Ramage, L., Nuki, G., and Salter, D. M. (2009). Signalling Cascades in Mechanotransduction: Cell-Matrix Interactions and Mechanical Loading. Scand. J. Med. Sci. Sports 19, 457–469. doi:10.1111/j.1600-0838.2009.00912.x
Ritchie, W. (2017). microRNA Target Prediction. Methods Mol. Biol. 1513, 193–200. doi:10.1007/978-1-4939-6539-7_13
Sekar, D. (2020). Implications of microRNA 21 and its Involvement in the Treatment of Different Type of Arthritis. Mol. Cel Biochem 476, 941. doi:10.1007/s11010-020-03960-y
Sharma, A., Jagga, S., Lee, S.-S., and Nam, J.-S. (2013). Interplay between Cartilage and Subchondral Bone Contributing to Pathogenesis of Osteoarthritis. Int. J. Mol. Sci. 14, 19805–19830. doi:10.3390/ijms141019805
Statello, L., Guo, C.-J., Chen, L.-L., and Huarte, M. (2021). Gene Regulation by Long Non-coding RNAs and its Biological Functions. Nat. Rev. Mol. Cel Biol. 22, 96–118. doi:10.1038/s41580-020-00315-9
Subramanian, A., Tamayo, P., Mootha, V. K., Mukherjee, S., Ebert, B. L., Gillette, M. A., et al. (2005). Gene Set Enrichment Analysis: a Knowledge-Based Approach for Interpreting Genome-wide Expression Profiles. Proc. Natl. Acad. Sci. 102, 15545–15550. doi:10.1073/pnas.0506580102
Sun, L., Luo, H., Bu, D., Zhao, G., Yu, K., Zhang, C., et al. (2013). Utilizing Sequence Intrinsic Composition to Classify Protein-Coding and Long Non-coding Transcripts. Nucleic Acids Res. 41, e166. doi:10.1093/nar/gkt646
Sun, H. B. (2010). Mechanical Loading, Cartilage Degradation, and Arthritis. Ann. N. Y Acad. Sci. 1211, 37–50. doi:10.1111/j.1749-6632.2010.05808.x
Tafer, H., and Hofacker, I. L. (2008). RNAplex: a Fast Tool for RNA-RNA Interaction Search. Bioinformatics 24, 2657–2663. doi:10.1093/bioinformatics/btn193
Tu, J., Huang, W., Zhang, W., Mei, J., and Zhu, C. (2020). The Emerging Role of lncRNAs in Chondrocytes from Osteoarthritis Patients. Biomed. Pharmacother. 131, 110642. doi:10.1016/j.biopha.2020.110642
van Hoolwerff, M., Metselaar, P. I., Tuerlings, M., Suchiman, H. E. D., Lakenberg, N., Ramos, Y. F. M., et al. (2020). Elucidating Epigenetic Regulation by Identifying Functional Cis-Acting Long Noncoding RNAs and Their Targets in Osteoarthritic Articular Cartilage. Arthritis Rheumatol. 72, 1845. doi:10.1002/art.41396
Vatner, S. F., Pachon, R. E., and Vatner, D. E. (2015). Inhibition of Adenylyl Cyclase Type 5 Increases Longevity and Healthful Aging through Oxidative Stress protection. Oxid. Med. Cel Longev. 2015, 250310. doi:10.1155/2015/250310
Villegas, V., and Zaphiropoulos, P. (2015). Neighboring Gene Regulation by Antisense Long Non-coding RNAs. Int. J. Mol. Sci. 16, 3251–3266. doi:10.3390/ijms16023251
Xiang, W., Zhang, J., Wang, R., Wang, L., Wang, S., Wu, Y., et al. (2018). Role of IFT88 in Icariin Regulated Maintenance of the Chondrocyte Phenotype. Mol. Med. Rep. 17, 4999–5006. doi:10.3892/mmr.2018.8486
Xu, T., Yang, K., You, H., Chen, A., Wang, J., Xu, K., et al. (2013). Regulation of PTHrP Expression by Cyclic Mechanical Strain in Postnatal Growth Plate Chondrocytes. Bone 56, 304–311. doi:10.1016/j.bone.2013.06.027
Xu, X.-W., Zhou, X.-H., Wang, R.-R., Peng, W.-L., An, Y., and Chen, L.-L. (2016). Functional Analysis of Long Intergenic Non-coding RNAs in Phosphate-Starved rice Using Competing Endogenous RNA Network. Sci. Rep. 6, 20715. doi:10.1038/srep20715
Xu, B., Xing, R., Huang, Z., Yin, S., Li, X., Zhang, L., et al. (2019). Excessive Mechanical Stress Induces Chondrocyte Apoptosis through TRPV4 in an Anterior Cruciate Ligament-Transected Rat Osteoarthritis Model. Life Sci. 228, 158–166. doi:10.1016/j.lfs.2019.05.003
Yan, P., Luo, S., Lu, J. Y., and Shen, X. (2017). Cis- and Trans-acting lncRNAs in Pluripotency and Reprogramming. Curr. Opin. Genet. Develop. 46, 170–178. doi:10.1016/j.gde.2017.07.009
Yang, K., Wu, Y., Cheng, P., Zhang, J., Yang, C., Pi, B., et al. (2016). YAP and ERK Mediated Mechanical Strain-Induced Cell Cycle Progression through RhoA and Cytoskeletal Dynamics in Rat Growth Plate Chondrocytes. J. Orthop. Res. 34, 1121–1129. doi:10.1002/jor.23138
Yao, X., Zhang, J., Jing, X., Ye, Y., Guo, J., Sun, K., et al. (2019). Fibroblast Growth Factor 18 Exerts Anti-osteoarthritic Effects through PI3K-AKT Signaling and Mitochondrial Fusion and Fission. Pharmacol. Res. 139, 314–324. doi:10.1016/j.phrs.2018.09.026
Yao, X., Sun, K., Yu, S., Luo, J., Guo, J., Lin, J., et al. (2021). Chondrocyte Ferroptosis Contribute to the Progression of Osteoarthritis. J. Orthop. Translat. 27, 33–43. doi:10.1016/j.jot.2020.09.006
Yu, G., Wang, L.-G., Han, Y., and He, Q.-Y. (2012). clusterProfiler: an R Package for Comparing Biological Themes Among Gene Clusters. OMICS: A J. Integr. Biol. 16, 284–287. doi:10.1089/omi.2011.0118
Zhang J, J., Yin, M., Peng, G., and Zhao, Y. (2018). CRNDE: An Important Oncogenic Long Non-coding RNA in Human Cancers. Cell Prolif. 51, e12440. doi:10.1111/cpr.12440
Zhang L, L., Zhang, P., Sun, X., Zhou, L., and Zhao, J. (2018). Long Non-coding RNA DANCR Regulates Proliferation and Apoptosis of Chondrocytes in Osteoarthritis via miR-216a-5p-JAK2-STAT3 axis. Biosci. Rep. 38, BSR20181228. doi:10.1042/BSR20181228
Zhang, J., Hao, X., Chi, R., Qi, J., and Xu, T. (2021). Moderate Mechanical Stress Suppresses the IL-1beta-induced Chondrocyte Apoptosis by Regulating Mitochondrial Dynamics. J. Cel Physiol. 236, 7504. doi:10.1002/jcp.30386
Zhao, H., and Gong, N. (2019). miR-20a Regulates Inflammatory in Osteoarthritis by Targeting the IκBβ and Regulates NK-κB Signaling Pathway Activation. Biochem. Biophys. Res. Commun. 518, 632–637. doi:10.1016/j.bbrc.2019.08.109
Keywords: osteoarthritis, chondrocyte biology, mechanical stress, transcriptome, non-coding RNA
Citation: Zhang J, Hao X, Chi R, Liu J, Shang X, Deng X, Qi J and Xu T (2022) Whole Transcriptome Mapping Identifies an Immune- and Metabolism-Related Non-coding RNA Landscape Remodeled by Mechanical Stress in IL-1β-Induced Rat OA-like Chondrocytes. Front. Genet. 13:821508. doi: 10.3389/fgene.2022.821508
Received: 24 November 2021; Accepted: 17 January 2022;
Published: 03 March 2022.
Edited by:
Chi-Ming Wong, Hong Kong Polytechnic University, Hong Kong SAR, ChinaReviewed by:
Agnieszka Fiszer, Institute of Bioorganic Chemistry (PAS), PolandStephen J Bush, University of Oxford, United Kingdom
Xueqiang Wang, Shanghai University of Sport, China
Copyright © 2022 Zhang, Hao, Chi, Liu, Shang, Deng, Qi and Xu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Jun Qi, cWlqdW5AdGpoLnRqbXUuZWR1LmNu; Tao Xu, eHV0ZG9jQDE2My5jb20=
†These authors have contributed equally to this work