Analysis of the lncRNA-Associated Competing Endogenous RNA (ceRNA) Network for Tendinopathy

Background We aimed to construct the lncRNA-associated competing endogenous RNA (ceRNA) network and distinguish feature lncRNAs associated with tendinopathy. Methods We downloaded the gene profile of GSE26051 from the Gene Expression Omnibus (GEO), including 23 normal samples and 23 diseased tendons. Differentially expressed mRNAs (DEmRNAs) and differentially expressed lncRNAs (DElncRNAs) were identified, and functional and pathway enrichment analyses were performed. Protein-protein interaction (PPI) network was constructed and further analyzed by module mining. Moreover, a ceRNA regulatory network was constructed based on the identified lncRNA–mRNA coexpression relationship pairs and miRNA–mRNA regulation pairs. Results We identified 1117 DEmRNAs and 57 DElncRNAs from the GEO data. The downregulated DEmRNAs were particularly associated with muscle contraction and muscle filament, while the upregulated ones were linked to extracellular matrix organization and cell adhesion. From the PPI network, 11 modules were extracted. Genes in MCODE 2 (such as TPM4) were significantly involved in cardiomyopathy, and genes in MCODE 4 (such as COL4A3 and COL4A4) were involved in focal adhesion, ECM-receptor interaction, and PI3K-Akt signaling pathway. The ceRNA network contained 7 lncRNAs (MIR133A1HG, LINC01405, PRKCQ-AS1, C10orf71-AS1, MBNL1-AS1, HOTAIRM1, and DNM3OS), 63 mRNAs, and 41 miRNAs. Downregulated lncRNA MIR133A1HG could competitively bind with hsa-miR-659-3p and hsa-miR-218-1-3p to regulate the TPM3. Meanwhile, MIR133A1HG could competitively bind with hsa-miR-1179 to regulate the COL4A3. Downregulated C10orf71-AS1 could competitively bind with hsa-miR-130a-5p to regulate the COL4A4. Conclusions Seven important lncRNAs, particularly MIR133A1HG and C10orf71-AS1, were found associated with tendinopathy according to the lncRNA-associated ceRNA network.


Introduction
Tendinopathy is characterized by a disorganized, haphazard healing response with no histological signs of inflammation. [1] Pathologies of the Achilles, patellar, and rotator cuff tendons are still common problems for recreational and competitive athletes and for individuals lacking any history of promoted physical activity. [2] A particular problem is rotator cuff tendinopathy (RCT), which is a common musculoskeletal disorder, accounting for 85% of shoulder joint pain. [3] Extrinsic factors result in a large number of acute tendon injuries, such as tendinopathy. Although in overuse syndromes, the combination of intrinsic factors may be the cause of tendinopathy, such as age-related changes in cell activity and external factors.
Recently, research focus has been placed on long noncoding RNAs (lncRNAs), which are a category of RNA molecules of over 200 nucleotides in length with an inconspicuous open reading frame. [4] Previous studies have reported that many human diseases, such as orthopedic and cancer diseases, were associated with deficiency, mutation, or overexpression of lncRNAs. Salmena et al. [5] proposed the competing endogenous RNA (ceRNA) hypothesis in 2011. Subsequently, it was asserted that ceRNAs play roles in orthopedic diseases, such as RCT and osteoarthropathy (OA). Ge et al. [3] identified several lncRNAs and mRNAs, such as COL11A1 and COL2A1, which show differential expression in RCT, suggesting their relation to the process of this disorder. In addition, Shen et al. [6] reported that SNHG5 refers to the mechanism of OA via acting as a ceRNA to competitively sponge miR-26a, thereby regulating the expression of SOX2. Moreover, Li et al. [7] demonstrated that the lncRNA XIST could increase OA chondrocytes proliferation and increase apoptosis via the miR-211/CXCR4 axis. erefore, the lncRNA XIST is a valuable curative target for treating OA. ese researches express that other lncRNAs might also occupy vital roles in the etiopathogenesis of orthopedic diseases.
In the study of Jelinsky et al. [8], the differentially expressed mRNAs (DEmRNAs) between normal and tendinopathy samples (mainly RCT tear and extensor carpi radialis brevis tear) were revealed, and pathway analysis was conducted, and the results showed that tendinopathy samples lacked appreciable inflammatory response and had altered expressions of extracellular matrix proteins. ereafter, Cai et al. downloaded the GSE26051 that was deposited by Jelinsky et al., the potential regulatory microRNA of DEmRNAs was further predicted, and especially, miR-499 was found to regulate specific genes, including CUGBP2 and MYB [9]. Previous studies contribute to the understanding of molecular mechanisms of tendinopathy. However, to the best of our knowledge, the roles of ceRNAs in the pathologies of tendinopathy have seldomly been investigated. In order to further elucidate the roles of lncRNAs in tendinopathy, we downloaded GSE26051 that deposited by Jelinsky et al. from Gene Expression Omnibus (GEO). e differentially expressed mRNAs (DEmRNAs) and differentially expressed lncRNAs (DElncRNAs) in tendinopathy samples were identified, followed by functional and pathway enrichment analyses. Protein-protein interaction (PPI) was constructed for module analysis. Moreover, an lncRNAassociated ceRNA regulatory network was constructed based on the lncRNA-mRNA coexpression relationship pairs and miRNA-mRNA regulation pairs. Hopefully, the results of our study could assist in understanding the roles of lncRNAs in tendinopathy and provide valuable treatment targets.

Affymetrix Microarray Dataset.
Gene expression profilings of tendinopathy were retrieved from National Center for Biotechnology Information (NCBI) GEO (https://www. ncbi.nlm.nih.gov/geo/) by searching the keywords of "Homo sapiens" and "chronic tendon injuries". e inclusion criteria of datasets were as follows: samples were collected from patients with chronic tendon injuries and contained normal samples; the number of overall samples was not less than 20. One dataset GSE26051 [8] meeting the above criteria was downloaded from GEO [10] (http://www.ncbi.nlm.nih.gov/ geo/), which was based on the GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array platform.
is dataset consisted of 46 specimens including 23 biopsies of diseased tendons and 23 biopsies of normal-appearing tendons from 23 patients who were undergoing surgical procedures for the treatment of chronic tendinopathy.

Data Preprocessing and Annotation.
We analyzed the derived genetic data using the affy package in R [11] (version 1.50.0, http://www.bioconductor.org/packages/release/bioc/ html/affy.html). en, we used the RMA (robust multiarray average) method to perform background correction, data normalization, and expression calculation. We obtained the expression values of probes in the normal and diseased tendons.
Subsequently, we obtained the sequence corresponding to the HG-U133_Plus_2 probe from the platform annotation file. We downloaded the current human reference genome (GRCh38) from the GENCODE database [12] (https://www.gencodegenes.org/releases/current.html).
en, we used SeqMap [13] to map all of the probe sequences to the reference genome. e unique mapping probes were retained, and then, the gene corresponding to each probe was obtained on the basis of the human gene annotation file (Release 25) provided by GENCODE according to the genomic features such as position in chromosome, sense, or antisense direction. Protein-coding probes were considered as probes of mRNAs, and probes with the annotation features of antisense, sense_intronic, lincRNA, sense_overlapping, or processed_transcript were considered as probes of lncRNAs.

GO BPs and Pathway Enrichment Analyses for
DEmRNAs. Analysis of the GO BP [16] and KEGG pathway enrichment annotations was performed for DEmRNAs on basis of e Database for Annotation, Visualization, and Integrated Discovery [12] (DAVID, version 6.8, https://david. ncifcrf.gov/) [17,18], with p < 0.05 and count ≥2 as the thresholds.
Moreover, GO and KEGG pathway enrichment analyses were carried out using the online tool Metascape (http:// metascape.org). [19] Parameters were set to min overlap � 3; p value � 0.01; min enrichment � 1.5; and enrichment factor >1.5. e enrichment factor is the ratio between the observed count and the accidentally expected count. After obtaining the terms that fulfill the above criteria, clustering analysis was carried out according to the similarity of genes (similarity >0.3) enriched in each term. e most statistically significant term with the minimum p value was considered as a representative of that cluster. e top 20 most significant clusters were retained for bar graph display. e online tool Metascape was applied to mine the PPI relationships of the above-mentioned differentially expressed genes. en, the databases BioGRID [20], inweb_ IM [21], and OmniPath [22] were used for protein-protein interaction enrichment analysis, using the default parameters (min network size � 3; max network size � 500). After getting the PPI pairs, the network was visualized using the software of Cytoscape (version 2.1.6, http://apps.cytoscape.org/), and the plugin CytoNCA (http://apps.cytoscape.org/apps/cytonca) [23] was used to analyze the degree of network nodes. e parameter was set without weighting.
rough ranking the connectivity of each node, we obtained the important nodes in the PPI network, namely the hub proteins. Furthermore, we used the Metascape tool to mine PPI network modules based on the molecular complex detection (MCODE) algorithm. [24] At the same time, we performed GO BP (GO biological process) and KEGG pathway enrichment analyses for each module.

Analysis of Coexpression of lncRNA and mRNA and Pathway Enrichment Analysis of lncRNA.
e Pearson correlation coefficients of the above-mentioned DEmRNAs and DElncRNAs were calculated, and the correlation test was carried out using the matched mRNA and lncRNA data. In addition, Benjamini-Hochberg (BH) procedure was used to obtain the adjusted p value [25]. To construct the subsequent ceRNA network, we focused on the pairs with correlation coefficient r >0.7 and adjusted p value < 0.05. e KEGG pathway enrichment analysis was performed for the target genes of each lncRNA that was enriched by using clusterProfiler in R (version: 3.8.1, http://bioconductor.org/ packages/release/bioc/html/clusterProfiler.html) [26], and the adjusted p value was obtained by BH procedure. Significantly enriched pathways were obtained with the threshold of the adjusted p value < 0.05. e top ten lncRNA pathways were displayed in a bubble chart.

Construction of the ceRNA Network.
Based on the miRNA-mRNA and lncRNA-miRNA relationship pairs obtained as described above, we first selected the lncRNA-miRNA-mRNA relationship pairs regulated by the same miRNA and then performed further lncRNA-miRNA-mRNA screening by combining the positive coexpression relationship between mRNA and lncRNA (correlation coefficient r > 0.7). e lncRNAs and mRNAs regulated by the same miRNAs in the ceRNA network that have a positive coexpression relationship were considered as ceRNAs. Finally, we used the Cytoscape plugin CycloNCA to analyze the node degree, and the parameter was set without weighting. A higher connection degree reflects the greater importance of the node in the network.

Identification of DEmRNAs and DElncRNAs in Tendinopathy.
A total of 37,666 probes were obtained through data preprocessing and standardization, which is shown in Appendix Table 1. e box diagram of the gene expression distribution of each sample is shown in Figure 1(a), and the median value of sample gene expression was basically at the same level, which could be directly used for subsequent difference analysis. When compared with normal controls, a total of 1507 probes were regarded as differentially expressed in tendinopathy samples, which corresponded to 1117 DEmRNAs (including 702 upregulated DEmRNAs and 415 downregulated ones) and 57 DElncRNAs (including 29 upregulated DElncRNAs and 28 downregulated ones) ( Table 1). Twenty-seven lncRNAs were screened from the featured lncRNAs in the LncBook database or lncRNAs recorded in the LNCiperia database (Appendix Table 2). Based on the expressions of the DEmRNAs and 27 functional lncRNAs, the heatmaps of DEmRNAs (Figure 1(b)) and functional lncRNAs (Figure 1(c)) showed that the samples were clustered in two different directions.

Gene Function and Pathway Enrichment Analyses.
Totally, 148 GO BPs and 16 KEGG pathways were highly enriched with a cutoff of p value < 0.05 for the upregulated DEmRNAs, while 79 GO BPs and 19 KEGG pathways were obtained for the downregulated DEmRNAs. e top 10 GO BPs and KEGG pathways (ranked by p value) are shown in Figure 2. Figure 2(a) reveals that the downregulated DEmRNAs were particularly associated with muscle contraction and muscle filament, while the upregulated DEmR-NAs were particularly associated with extracellular matrix organization and cell adhesion. Furthermore, the KEGG pathway enrichment analysis findings indicated that the most remarkable pathway was focal adhesion (Figure 2(b)).
Besides, through the KEGG pathways enrichment analysis and GO BPs in the Metascape database, five pathways involved in the downregulated DEmRNAs included muscle structure development, actin filament-based movement, muscle system process, actin filament-based process, and skeletal muscle contraction. Besides, four pathways shown to be related to the upregulated DEmRNAs were ossification, blood vessel development, extracellular matrix organization, and skeletal system development (Figure 3(a)).
Moreover, term-term interaction network (Figures 3(b)) showed that the upregulated and downregulated DEmRNAs participate in different functions. e network consists of 172 nodes, representing 172 terms, and 947 edges were enriched, which represent the similarity score between terms.
Metascape obtained the similarity score between the two terms. e larger the score, the thicker the connection lines.

Construction of PPI Network and Module Mining
Analysis. A total of 2095 PPI pairs, including 440 DEmR-NAs, were identified to construct the PPI network ( Figure 4(a)). In the PPI network, the top 10 genes ranked by node degree were ACTB, CDK1, ACTN2, ACTA1, TPM1, ACTA2, ACTN4, FLNA, ENO1, and ACTN1 in order. Moreover, 11 module substructures were distinguished in the PPI network (Figures 4(b)). e pathway enrichment analysis for the genes in the 11 Table 3).

Coexpression Analysis of DEmRNAs and DElncRNAs.
A total of 1851 significant coexpression relationships were screened with a cutoff of p < 0.05, which included 422 mRNAs and 22 lncRNAs. KEGG pathways were enriched for the target genes of 11 lncRNAs among the 22 lncRNAs, and the top ten KEGG pathways are shown in Figure 5. e results showed that the target genes of C10orf71−AS1, LINC00844, LINC01405, MIR133A1HG, PCAT7, and PRKCQ−AS1 were particularly associated with cardiac muscle contraction, followed by adrenergic signaling in cardiomyocytes and calcium signaling pathway ( Figure 5).

Discussion
In our study, our purpose was to distinguish potential lncRNAs and mRNAs involved in tendinopathy. Overall, data from GEO revealed 1117 DEmRNAs and 57    Note. e number outside the bracket is the horizontal number of genes, and the number inside the bracket is the number of corresponding probes.

Genetics Research 5
DElncRNAs. Next, we identified the function and signaling pathways, where these DEmRNAs and lncRNAs were embroiled. What's more, we established a ceRNA coexpression network that contained 7 lncRNAs, 63 mRNAs, and 41 miRNAs using bioinformatic tools.
en, Ge et al. [3] constructed a coexpressed network comprising the vital genes including NONHSAT209114.1, ENST00000577806, NONHSAT168 464.1, PLK2, TMEM214, and IGF2 for tendinopathy. Furthermore, a mechanism study demonstrated that H19 promotes tenogenic differentiation both in vitro and in vivo by targeting miR-29b-3p and activating TGF-β1 signaling [30]. As previous study indicated that lncRNAs also played a key role in male infertility [31], Rotondo et al. [32] reported that lncRNA H19 was dysregulated in male infertility by epigenetic impairments. e defective marks of H19 could be potentially employed as useful tools in clinical practice for assessing male infertility [33]. Nafiseh et al [34] found that SLC7A11-AS1 played a potential role in the pathophysiology of male infertility associated with varicocele. ese reports indicated that lncRNA expression could regulate the development of some diseases. Further experimental and functional studies need to be performed for exploring the role of lncRNAs. When such studies are done, these findings should help guide them in the future. e important BPs and Genetics Research pathways in tendinopathy were identified using GO and KEGG pathway analyses on basis of 1117 DEmRNAs. Most of the BPs and pathways play a crucial role in tendinopathy, for example, muscle contraction, muscle filament, matrix organization, cell adhesion, and focal adhesion. [35] Besides, eleven intersecting lncRNAs that overlap between the GO and KEGG pathway analyses were identified in tendinopathy. e results showed that C10orf71−AS1, LINC00844, LINC01405, MIR133A1HG, PCAT7, and PRKCQ−AS1 were particularly associated with cardiac muscle contraction, followed by adrenergic signaling in cardiomyocytes. Previous studies indicated that cardiac muscle contraction is involved in diabetes and ischemia/reperfusion oxidative injury. [36] Moreover, a total of 2095 PPI pairs, including 440 DEmRNAs, were contained in the PPI network. In total, 10 hub genes, namely, ACTB, CDK1, ACTN2, ACTA1, TPM1, ACTA2, ACTN4, FLNA, ENO1, and ACTN1, were screened from the PPI network. Tsai et al. [37] indicated that geranylgeranyl pyrophosphate could prevent the adverse effect of simvastatin in tendon cells through downregulating CDK1 expression. In contrast, these other hub genes have not been reported to be related to tendinopathy. Moreover, 11 modules were extracted from the PPI network. e results showed that genes (DES, MYH6, MYL2, MYL3, TPM1, TPM3, TPM4, and TTN) in MCODE 2 were significantly involved in dilated cardiomyopathy and hypertrophic cardiomyopathy (HCM). It has been reported that rupture of the distal biceps tendon was found in 33.3% of patients with wild-type transthyretin amyloidosis cardiomyopathy [38]. Meanwhile, genes (COL4A1, COL4A2, COL4A3, COL4A4, COL6A1, COL6A2, CTNNB1, and PDGFB) in MCODE 4 were involved in focal adhesion, ECM-receptor interaction, and PI3K-Akt signaling pathway. e transcription factor scleraxis (Scx) plays a critical role in tenocyte mechanotransduction, and focal adhesions and extracellular matrixreceptor interaction were significantly downregulated in Scx knockdown tenocytes [39]. Inhibition of adipogenesis of tendon stem cells and fatty infiltration in injury tendon could promote biomechanical properties and decrease rupture risk of injury tendon by downregulating PTEN/ PI3K/AKT signaling [40].
us, genes in MODE 2 and MODE 4 might be associated with tendinopathy by regulating cardiomyopathy, focal adhesion, ECM-receptor interaction, and PI3K-Akt signaling pathway.
Moreover, on basis of the lncRNA-mRNA coexpression relationship pairs and miRNA-mRNA regulation pairs, a lncRNA-associated ceRNA network was constructed. In the ceRNA network, there were 7 lncRNAs, 63 mRNAs, and 41 miRNAs. Downregulated lncRNA MIR133A1HG could competitively bind with hsa-miR-659-3p and hsa-miR-218-1-3p to regulate the expression of TPM3 and might influence      the cardiomyopathy. Meanwhile, MIR133A1HG could competitively bind with hsa-miR-1179 to regulate the expression of COL4A3, and downregulated C10orf71-AS1 could competitively bind with hsa-miR-130a-5p to regulate the expression of COL4A4, thereby might regulate the focal adhesion, ECM-receptor interaction, and PI3K-Akt signaling pathway. However, our study still has some limitations. Firstly, the main limitation of the study is the lack of validation in vitro/ vivo as well as functional validation on the role of lncRNAs and mRNAs in tendinopathy with cell lines. en, the sample size of our study was small in the computational analyses, and bigger samples were necessary to validate these lncRNAs' functions in tendinopathy.
Data Availability e dataset used and/or analyzed during the current study are available from the corresponding author upon reasonable request.

Ethical Approval
Ethical approval was not required for this study.

Conflicts of Interest
e authors declare that there are no conflicts of interest.

Authors' Contributions
Xiao-Ning Liu and Bing-Zhe Huang were responsible for the conception and design of the research and drafting of the manuscript. Jing-Jing Yang and Xiao-Ming Dong performed the data acquisition. Jing-Jing Yang and Zuan Zhong performed the data analysis and interpretation. Bing-Zhe Huang participated in the design of the study and performed the statistical analysis. All authors have read and approved the manuscript. Acknowledgments e study was supported by the Special Foundation for Science and Technology Innovation of Jilin Province (NO. 20200201566JC). Table 1 shows data preprocessing and standardization. Table 2 shows featured lncRNAs in the LncBook database or lncRNAs recorded in the LNCiperia database. Table 3 shows pathway enrichment analysis for the genes in the 11 modules. Supplementary Materials (Supplementary  Materials)